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Abstract 

We provide a comprehensive survey of splitting and composition meth- 
ods for the numerical integration of ordinary differential equations (ODEs). 
Splitting methods constitute an appropriate choice when the vector field 
associated with the ODE can be decomposed into several pieces and each 
of them is integrable. This class of integrators are explicit, simple to im- 
plement and preserve structural properties of the system. In consequence, 
they are specially useful in geometric numerical integration. In addition, 
the numerical solution obtained by splitting schemes can be seen as the 
exact solution to a perturbed system of ODEs possessing the same geomet- 
ric properties as the original system. This backward error interpretation 
has direct implications for the qualitative behavior of the numerical so- 
lution as well as for the error propagation along time. Closely connected 
with splitting integrators are composition methods. We analyze the order 
conditions required by a method to achieve a given order and summarize 
the different families of schemes one can find in the literature. Finally, 
we illustrate the main features of splitting and composition methods on 
several numerical examples arising from applications. 
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1 Introduction by examples 

The basic idea of splitting methods for the time integration of ordinary differ- 
ential equations (ODEs) can be formulated as follows. Given the initial value 
problem 

x' = f(x), x = x{0) e M D (1) 
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with / : R D — > R D and solution (p t (x ), let us suppose that / can be expressed 
as / = 5^2=1 / f° r certain functions /W : M' — ► M. D , in such a way that the 
equations 

x' = f®(x), x = x(0)eR D , i = l,...,m (2) 

can be integrated exactly, with solutions x(h) = ip$(xo) at t = h, the time step. 
Then, by combining these solutions as 

[ml [2] [1] /q\ 

Xh = pi 1 ° ■ ■ ■ ° v>h ° <p l h ( 3 ) 

and expanding \ m t° Taylor series, one finds that Xh( x o) = (Phfao) + C(^ 2 )> so 
that Xh provides a first-order approximation to the exact solution. As we will 
see, it is possible to get higher order approximations by introducing more maps 
with additional coefficients, ^.. h , in the previous composition ([3]). 

One thus may say that splitting methods involve three steps: (i) choosing 
the set of functions such that / = ( n ) solving either exactly or 

approximately each equation x' = /^(x); and (iii) combining these solutions 
to construct an approximation for ([I]). One obvious requirement is that the 
equations x' = f^(x) should be simpler to integrate than the original system. 

Some of the advantages that splitting methods possess can be summarized 
as follows: 

• They are usually simple to implement. 

• They are, in general, explicit. 

• Their storage requirements are quite modest. The algorithms are sequen- 
tial and the solutions at intermediate stages are stored in the solution 
vectors. This property can be of great interest when they are applied to 
partial differential equations (PDEs) previously semidiscretized. 

• There exist in the literature a large number of specific methods tailored 
for different structures. 

• They preserve structural properties of the exact solution, thus conferring 
to the numerical scheme a qualitative superiority with respect to other 
standard integrators, especially when long time intervals are considered. 
Examples of these structural features are symplecticity, volume preser- 
vation, time-symmetry and conservation of first integrals. In this sense, 
splitting methods constitute an important class of geometric numerical 
integrators. 

Let us give more details on this last item. Traditionally, the goal of numer- 
ical integration of ODEs consists in computing the solution to the initial value 
problem ([T|) at time tiy = Nh with a global error \\xj\f — x(tjv)|| smaller than a 
prescribed tolerance and as efficiently as possible. To do that one chooses the 
class of method (one-step, multistep, extrapolation, etc.), the order (fixed or 
adaptive) and the time step (constant or variable). In contrast, with a geomet- 
ric numerical integrator one typically fix a (not necessarily small) time step and 
compute solutions for very long times for several initial conditions, in order to 
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get an approximate phase portrait of the system. It turns out that, although 
the global error of each trajectory may be large, the phase portrait is in some 
sense close to that of the original system. 

The aim of geometric numerical integration is thus to reproduce the qual- 
itative features of the solution of the differential equation which is being dis- 
cretised, in particular its geometric properties \17\ I41j . The motivation for 
developing such structure-preserving algorithms arises independently in areas 
of research as diverse as celestial mechanics, molecular dynamics, control the- 
ory, particle accelerators physics, and numerical analysis [HI [44] [57], [58l US] . 
Although diverse, the systems appearing in these areas have one important 
common feature. They all preserve some underlying geometric structure which 
influences the qualitative nature of the phenomena they produce. In the field 
of geometric numerical integration these properties are built into the numerical 
method, which gives the method an improved qualitative behaviour, but also 
allows for a significantly more accurate long-time integration than with general- 
purpose methods. In addition to the construction of new numerical algorithms, 
an important aspect of geometric integration is the explanation of the relation- 
ship between preservation of the geometric properties of the scheme and the 
observed favorable error propagation in long-time integration. 

Before proceeding further, let us introduce at this point some splitting meth- 
ods and illustrate them on simple examples. 



Example 1: Symplectic Euler and leapfrog schemes. Suppose we have 
a Hamiltonian system of the form H(q,p) = T(p) + V(q), where g £ l d are 
the canonical coordinates, p £ M rf are the conjugate momenta, T represents the 
kinetic energy and V is the potential energy. Then the equations of motion 
read [36] 

q' = T p (p), p' = -V q (q), (4) 

where T p and V q denote the vectors of partial derivatives. Equations (pi]) can be 
formulated as P with x = (q,p) T , f(x) = (T p , -V q ) T = JVH(x) and D = 2d. 
Here J denotes the 2d x 2d canonical symplectic matrix 



J 



I d 
-Id o 



and Id stands for the d-dimensional identity matrix. In this case the exact flow 
(ft is symplectic [TJ. The simple Euler method applied to this system provides 
the following first order approximation for a time step h: 

q n+ i = q n + hT p {p n ) , . 

p n+ i = p n -hV q (q n ). 

On the other hand, if we consider H as the sum of two Hamiltonians, the first 
one depending only on p and the second only on q, the corresponding Hamilton 
equations 

q ' = T ^ P) and q ' = ° (6) 

p' = and p> = -V q (q) (6) 
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with initial condition (qo,Po) can be readily solved to yield 



[T] : q(f) = Q0 + tT p ( P0 ) [y] q(t) = q 

^ P(t) = po ^ P(t) = P0-tV q (q ), 

(7) 

respectively. Composing the time t = h flow ip^ (from initial condition 

[71 

{Qn,Pn)) followed by (p h , gives the scheme 

Xh = lfh 0lfh ■ q n +i = q n + hT p (p n+1 ). (8) 

Since it is a composition of the flows of two Hamiltonian systems and in addition 
the composition of symplectic maps is again symplectic, Xh is a symplectic 
integrator, which can be called appropriately the symplectic Euler method. It 

[VI [71 

is of course also possible to compose the maps in the opposite order, (p l h o (p l h , 
thus obtaining another first order symplectic Euler scheme: 

Xk = lPh 0(Ph ■ Pn + 1 = Pn-hV q (q n+1 ). ^ 

One says that ([9]) is the adjoint of Xh- Yet another possibility consists in using 
a 'symmetric' version 

c [2] _ [V] [T] [V] / ln s 

which is known as the Strang splitting [77], the leapfrog or the Stormer-Verlet 



[2l 

method [85] , depending on the context where it is used. Observe that S h = 
Xh/2 ° Xh/2 anc ^ ^ ^ S a ^ so symplectic and second order. 

Example 2: Harmonic oscillator. Let us consider now the Hamiltonian 
function H(q,p) = + q 2 ), where now q,p € M. Then the corresponding 
equations (j3|) are linear and can be written as 



x 



' N r/ 1\ ( M ( q 

-(A + B)x. (11) 



/ + V -1 



V 



This system has periodic solutions for which the energy H is conserved. In 
addition, it is area preserving and time reversible. The numerical solution 
obtained by the Euler scheme (|5|) reads 



qn+i \ = ( 1 h \ / q n 

Pn+1 ) \ ~h I ) \ p n 

whereas the symplectic Euler method Q leads to 



(12) 



q n+ i \ I h \ q n \ =e kB eh A Qn\ (13) 

Pn+1 J \ ~h I - h Z J \ p n J \ Pn J 

Both render first order approximations to the exact solution, which can be 
expressed as x(t) = e h ^ A+B ^XQ, but there are important differences between 
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them. First, the map (|13j) is area preserving (because it is symplectic), in 
contrast with <\12\) . Second, the approximation obtained by the symplectic 
Euler scheme verifies 

^(p 2 n+ l + hp n+ iq n+ i + q 2 +1 ) = -{p 2 + hp n q n + q 2 n ). 

Third, it can be shown that (|13|) is the exact solution at t = h of the perturbed 
Hamiltonian system 

T ~w i \ 2 arcsin(/i/2) . 9 9 

H(q,p,h) = h ^ r L J (p 2 + hpq + q 2 ) (14) 

= \{P 2 + q 2 )+h {^pq+^h{p 2 +q 2 ) + ---^j . 

In other words, the numerical approximation (|13p . which is only of first order 
for the exact trajectories of the Hamiltonian H(q,p) = \{p 2 + q 2 ), is in fact the 
exact solution of the perturbed Hamiltonian (|14[) . 

How these features manifest in actual simulations? To illustrate this point 
we take initial conditions (qo,Po) = (4, 0) and integrate with a time step h = 0.1. 
Figure Q] shows the first five numerical approximations obtained by the Euler 
method (|12p and the symplectic Euler scheme (|13p in the left panel, and the 
results for the first 100 steps in the right panel. It is clear that for one time 
step there are not significant differences between the standard Euler and the 
symplectic Euler methods, but the picture is completely different for longer 
integrations, where the superiority of the splitting symplectic method is evident. 
Note that the numerical solution it provides evolves on a slightly perturbed 
ellipse. 



Example 3: The 2-body problem (Kepler problem). The motion of two 
bodies attracting each other through the gravitational law can be described by 

Qi=~ ( ^tW . *' = M (15) 

in conveniently normalized coordinates in the plane of motion. This system has 
a number of characteristic geometric properties. First, equations (fT5|) can be 
derived from the Hamiltonian function 



H(q,p) = T(p) + V(q) = ±(p 2 + p 2 ) - i r = y/qj+q 2 . 



(16) 



Second, it is invariant under continuous translations in time and rotations in 
space, and thus both the Hamiltonian and the angular momentum L = q\P2 — 
q2P\ are preserved. In addition, the so-called Laplace-Runge-Lenz vector is 
also constant along their solutions, due to the fact that the symmetry group of 
this problem is the group of four-dimensional real proper rotations SO(4) |36j . 
For the numerical integration of this problem we choose as initial value 



gi(0) = l-e, g 2 (0)=0, pi(0) = 0, P2 (0) = j\+±, (17) 
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h=1/10 



h=1/10 




Figure 1: Numerical integration of the harmonic oscillator using the Euler 
method (white circles) and the symplectic Euler method (black circles) with 
initial condition (qo,Po) = (4,0) and time step h = j^. The left panel shows 
the results for the first 5 steps, whereas the right panel shows the results for 
the first 100 steps. The exact solution corresponds to the solid line. 

where < e < 1 represents the eccentricity of the orbit. In this case the total 
energy is H = Hq = —1/2, the period of the solution is 2ir, the initial condition 
corresponds to the pericenter and the major semiaxis of the ellipse is 1. 

Figure [2] shows some numerical solutions obtained with schemes (jSJ) and (jSJ) 
for the initial conditions (|17|) with eccentricity e = 0.6. The left panel shows 
the results for the integration of 3 periods with time step h = j^. As in 
the previous example, the explicit Euler method provides an approximate orbit 
that spirals outwards, whereas the symplectic Euler scheme merely distorts the 
ellipse, but also exhibits a precession effect. To better illustrate this effect, 
we repeat the experiment for a longer interval (15 periods) and a larger time 
step (h = T^j) in the right panel. The explanation of these phenomena can be 
formulated as follows. On the one hand, the symplectic Euler method exactly 
conserves the angular momentum. On the other hand, the numerical solution 
it provides can be seen as the exact solution of a slightly perturbed Kepler 
problem, and thus SO(4) is no longer the symmetry group of the problem, so 
that the Laplace-Runge-Lenz vector is not preserved and the trajectories are 
not closed anymore. 

Next we check how the error in the preservation of energy and the global 
error in position propagates with time. For comparison, we also include the 
results obtained with a Runge-Kutta method of order 2 (Heun's method) and 
the leapfrog scheme (llOj) . We now consider e = 1/5 and integrate for 500 
periods. The step size is h = ^|^q in all cases, except for the Heun method 
which uses h = ^ instead. In this way all methods require the same number 
of force evaluations (Heun's method computes twice the force per step). The 
corresponding results are shown in Figure [3] in a logdog scale. Notice that the 
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Figure 2: Numerical integration of the 2-body problem using the Euler method 
(white circles) and the symplectic Euper scheme (black circles) for the initial 
conditions (|17p with eccentricity e = 0.6. The left panel shows the results for 
h = and the first 3 periods and the right panel shows the results for h = ^ 
and the first 15 periods. 

average error in energy does not grow for the split symplectic methods and 
the error in positions grows only linearly with time, in contrast with Euler and 
Heun schemes. The Stormer-Verlet integrator provides more accurate results 
than the Heun method with the same computational cost. 

A collection of (additional) examples. Splitting methods constitute an 
important tool in different areas of science. In addition to Hamiltonian systems, 
they can be successfully applied in the numerical study of Poisson systems, sys- 
tems possessing integrals of the motion (such as energy and angular momen- 
tum) and systems with (continuous, discrete and time-reversal) symmetries. As 
a matter of fact, splitting methods have been designed (often independently) 
and extensively used in fields as distant as molecular dynamics, simulation of 
storage rings in particle accelerators, celestial mechanics and quantum physics 
simulations. To see why this is so, next we collect a number of differential 
equations which appear in different contexts ranging from Celestial Mechanics 
to electromagnetism and Quantum Mechanics. These examples also try to illus- 
trate the fact that very often one particular equation can be split into different 
ways and the most appropriate methods may depend on the split chosen. 
We (arbitrarily) classify our examples into three different categories. 

1. Hamiltonian systems. 

(a) Generalized harmonic oscillator (M,N £ JR rfx,i ): 

H = X - V T My + -q T Nq 



(18) 



ERROR ENERGY: e=1/5 



ERROR POSITION: e=1/5 




LOG(t) 




2.5 3 3.5 



Figure 3: Error growth in energy and position for the Kepler problem with 
e = 1/5 and 500 periods achieved by the first order symplectic Euler (EulerSI) 
and the second order Stormer-Verlet integrator (SI2). For comparison, we also 
include the first order Euler and the second order Heun (RK2) methods. The 
time step is adjusted in such a way that all methods use 1500 force evaluations 
per period. 



(b) Henon-Heiles Hamiltonian |43j : 



H = \{p\ +p\) + \{q\ + ql) + qfa ~ \ql- 



(19) 



(c) Perturbed Kepler problem. It models the dynamics of a satellite mov- 
ing into the gravitational field produced by a slightly oblate planet: 



(20) 



where e is typically a small parameter. When e = 0, the 2-body 
problem (fTBj) is recovered. 

(d) The gravitational iV-body problem (qt, pi £ M 3 , i = 1, . . . , N): 



H 



2 ^ 2m 

i=l 



T 
■Pi Pi 



N i-1 
i=2 j=l 



1j\ 



(21) 



(e) The motion of a charged particle in a constant magnetic field per- 
turbed by k electrostatic plane waves [2Tj : 



2 L 2 

2. More general dynamical systems. 
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H(q,p,t) = -p 2 + -q 2 + e^2cos(q - ant). 



(22) 



i=l 
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(a) The Volterra-Lotka problem |41| . 



d_ 

dt 



u 

V 



-2 
1 



u 

V 



+ 



uv u(v — 2) 

-uv v(l — u) 

(23) 

with first integral I(it, u) = log it — u + 2 log v — v. 
(b) The Lorenz system [52^ 139] (split into linear and non- linear parts): 



d 
dt 


X 




' -cr a " 


X 




" " 


y 




r -1 


y 


+ 


-x 


z 




b _ 


z 




Ox 


Here 


a, r, 


b > 


are constant. 


The values considered 



X 

y 

z 



(24) 



o- = 10, r = 28 and b = 8/3. 

(c) The ABC-flow (/ = /a + Jb + /cs but other splits are also possible 
(571): 



— (a;, y, z) = A(0, s'mx, cos x) + B(cos y, 0, siny) + C(sin z, cos z, 0). 

(25) 

3. Evolutionary PDEs. 

Although we are mainly concerned here with splitting methods applied to 
ODEs, it turns out that they can also be used in the numerical integration 
of certain partial differential equations. Specifically, a number of PDEs 
relevant in the applications, after an appropriate space discretization, 
lead to a system of ODEs which can be subsequently solved numerically 
by splitting methods. Among these equations, the following are worth to 
be mentioned. 



(a) The Schrodinger equation (h= 1): 

1 



2m 



V 2 + V(x) *(a?,t). 



(26) 



(b) The Gross-Pitaevskii equation [68] : 

1 



2m 



V 2 + V(x) + a\y(x,t)\ 2 )V(x,t) (27) 



(c) The Maxwell equations 



1 d 1 

— VxE, — E = -VxB, 
\i at e 



(28) 



where E(x,t), B(x,i) are the electric and magnetic field vectors, 
fi(x) is the the permeability and e(x) is the permittivity. 
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As we stated before, splitting methods form a subclass of geometric numer- 
ical integrators for various types of ODEs. The reason is easy to grasp from 
the examples analyzed before. Suppose that the flow of the original differen- 
tial equation (TjQ) forms a particular group of diffeomorphisms (in the case of 
Hamiltonian system, the group of symplectic maps). If / is conveniently split 
as / = Yli f (step (i) in the construction process of a splitting scheme) and 
the flows corresponding to each /M also belong to the same group of diffeomor- 
phisms in such a way that they can be explicitly obtained (step (ii)), then, by 
composing these flows (step (hi)) we get an approximation in the group, thus 
inheriting geometric properties of the exact solution. These considerations also 
hold (with some modifications) when the exact flow forms a semigroup or a 
symmetric space. 

With respect to steps (i) and (ii) before, several comments are in order. 
First, whereas for certain classes of ODEs, the splitting can be constructed sys- 
tematically for any /, in other cases no general procedure is known, and thus 
one has to proceed on a case by case basis. Second, sometimes a standard split- 
ting is possible for a certain /, but there exist other possible choices leading to 
more efficient schemes (we will see some examples in section |8]). Third, whereas 
the original system possesses several geometric properties which are interest- 
ing to preserve by the numerical scheme, different splittings preserve different 
properties and it is not always possible to find one splitting preserving all of 
them. These aspects have been analyzed in detail in [57], where a classification 
of ODEs and general guidelines to find suitable splittings in each case is pro- 
vided. Here, by contrast, we will concentrate on the third step of any splitting 
method: given a particular splitting, we will show how to combine the flows 
of the pieces /W to get efficient higher order approximations. In any case, the 
reader is referred to the excellent review paper [57] and the monographs [41^149]. 
where these and other issues, mainly in connection with geometric numerical 
integration, are thoroughly examined. 

2 Splitting and composition methods 

2.1 Increasing the order of an integrator by composition 

It is well known that numerical integrators of arbitrarily high order can be 
obtained by composition of a basic integrator of low order. Consider for instance 
the leapfrog scheme ()10p . which is a second-order integrator : M. 2d — > ~R 2d . 
Then, a 4th order integrator : M? d — ► M 2d can be obtained as 

4 4] =^«sgo«sS, with a = ^W' ^ = 1 - 2a - ( 29 ) 

More generally, if one recursively defines S^ 2k+2 ^ : M? d — > M. 2d for k = 1, 2, ... as 

~[2fc+2] _ <J2fc] q [2k] q [2k] /nM 

with 

a= 2- 2V(2fc+D ' /3 = 1 -2«> (31) 
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then the schemes 5^ 2fc ' are of order 2k (k > 1) |78^ I88| . We will prove later 
on this assertion. At this point it is useful to introduce the notion of adjoint 
of a given integrator iph [73]. By definition, this is the method ip^ such that 
i/jfr = ipZh- A method that is equal to its own adjoint is called self-adjoint or 
(time-) symmetric. In this case, ip-h ° iph = id. Since the leapfrog scheme (jlOp 
can be rewritten as 

sf=X h /2°xi/ 2 , (32) 

[21 

where Xh is the symplectic Euler method ([8]) , then S h is certainly time- 
symmetric. Actually, given any basic first order integrator \h ■ ^ D — ► ^ D 
for the ODE system ([Tjh the composition f|32|) is a time-symmetric method of 
order 2, and the other way around: any self-adjoint second order integrator 
can be expressed as (|32p . Furthermore, the integrators Si (k = 1,2,...) 
recursively defined by (I30p - (l3ip are time-symmetric methods of order 2k. In 
particular, if f(x) in the ODE (pQ) is split as 

m 

/(*)=£/ M o»o (33) 

i=l 

then, time-symmetric integrators <s! 2fe ' of order 2k can be constructed in this 
way by considering the basic first order integrator 

M [2] [1] (oa\ 

and its adjoint 

* -i [i] [2] H 
Xh = X- h = <P l h °V>h °"' 0( Ph ■ 

This general procedure of constructing geometric integrators of arbitrarily high 
order, although simple, presents some drawbacks. In particular, the resulting 
methods require a large number of evaluations and usually have large truncation 
errors. 

As we will show, efficient schemes can be built by considering more general 
composition integrators. First observe that the (2fc)th order integrators <?[ 2 ^ 
can be rewritten in the form 

il>h = Xa 2s h o X* 2s _ lh o • • • O X a 2 h ° X* ai h (35) 

with s = 3 fc_1 and some fixed coefficients (ai, . . . , «2s) G M. 2s . Then the idea is 
to consider composition integrators of the form (|35p with appropriately chosen 
coefficients (ai, . . . , a2s) G M. 2s . 

In the particular case where the ODE ([1]) is split in two parts / = + 

and Xh = ° > one can trivially check that the composition integrator (|35p 
can be rewritten as 

/ [b] [a] [b] [h] [a] [b] 

where &i = a\ and for j = 1, . . . , s, 

Oj = «2j-l + 02j, = "2j + 02J+1 (37) 
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(with a2s+i = 0). Conversely, any integrator of the form f|36|) satisfying that 

J2t=i a i = Ylt=l hi can be expressed in the form (f35j) with \h = <Jh ° • ^ or 
later reference, we state the following result, due to McLachlan |54j . 

Theorem 1 The integrator A36\) is of order r for ODEs of the form fl]) with 
f : M. D — ► M. D arbitrarily split as f = +f^ if and only if the integrator &35\) 
(with coefficients ctj obtained from jS7\ )) is of order r for arbitrary consistent 
integrators %h- 

2.2 Integrators and series of differential operators 

Before proceeding further with the analysis, let us relate generic numerical 
integrators with formal series of differential equations. This relationship will 
allow one to formulate in a rather simple way the conditions to be satisfied by 
an integration scheme to achieve a given order of consistency. 

First of all, let us recall that an integrator iph : M> D — ► M. D for the system 
([1]) is said to be of order r if for all x G W D 

^ h {x) = ip h {x) + 0{h r+l ) (38) 

as h — > 0, where cph is the h-Row of the ODE ([T]). 

It is well known that, for any smooth function g : M. D — ► R, it formally holds 
that [66] 

g(<p h (x)) = g(x) + V -.F n [g]{x) = exp(hF)[g](x), 
z — ' n! 

n>l 

where i 7 is the Lie derivative associated to the ODE system ([T|), i.e., the first 
order linear differential operator F acting on functions in C°°(M. D , M.) as follows. 
For each g G C°°(R D , R) and each x = (xi, . . . , x D ) G R D 




where /(x) = (/i(x), . . . , /d(x)) t . Motivated by this fact, we consider for a 
basic integrator x/i : — ^ the linear differential operators X n (n > 1) 
acting on smooth functions g G C°°(R Z? ,R) as follows: 

1 d n 

Xn[g]{x) = -— g{ Xh {x))\ h= ^ (40) 
so that formally g(xh( x )) = ^(^)b]( 2; )] where 

=/ + ^W„, (41) 

n>l 

and / denotes the identity operator. Thus, the integrator Xh is of order r if 

X n = -F n , 1 < n < r. 
nl 
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Alternatively, one can consider the series of vector fields 

(_-i\m+l 

Y(h) = £ h n Y n = ]og(x(h)) = — ( hx i + h2x 2 + ■ ■ • r ' . 

n>l m>l 

that is, 

Yn= Yl m Yl X ^ ' ' ' X ^ 

m>l jiH hjm=n 

so that X(h) = exp(Y(h)), and formally, g(xh(x)) = exp(Y(h))\g](x). Clearly, 
the basic integrator is of order r if 

Y\ = F, Y n = for 2 < n < r. 

For the adjoint integrator x% = XZ\i one obviously gets g(Xh( x )) — e~ Y (~ h ^ [g](x) 
This shows that Xh is time-symmetric if and only if Y(h) = hY\ + h s Y^ + • • • , 
and in particular, that time-symmetric methods are of even order. 

It is possible now to check that the symmetric integrators given by 

rol 

(|30|) - (f3"T|) are actually schemes of order 2k provided that S h is a symmetric 
second order integrator. Consider the series of differential operators 

F™(h) = hF + h^F^ + h 2k+3 F^% + ■■■ 
such that g(S [2k] (x)) = exp(F^ (h))[g](x) . Then one clearly has 

exp{F^ 2k+2 \h)) = exp{F^{ah))exp{F^(f3h))exp{F^(ah)) 
which implies 

F^ k+2 \h) = h (2a + P) F + h 2k+1 (2a 2k+1 + (3 2k+1 ) if^ + 0(h 2k+3 ), 

and thus 5^ fc+2 is of order 2k + 2 provided that Sf^ is of order 2k and a and (5 
satisfy the equations 

2a + [3 = I, 2a 2k+1 + (3 2k+1 = 0, 

whose unique real solution is given by (f3T1) . 

In the general case, for the composition method (|35p we have 

g(^ h (x)) = ^(h)[g\(x), 
where ^(h) = I + h^i + h 2 ^ 2 + ■ ■ ■ is a series of differential operators satisfying 
= X(-a l h)~ 1 X(a 2 h) ■ ■ ■ X(-a 2s -ih)- l X(a 2s h), (42) 
where the series X(h) is given by (|40p ~ (|4ip . and 

X(h)- 1 =I+Y1 ("l) m+1 {hXi + h 2 X 2 + •••)". (43) 

m>l 
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Thus, the order of a composition integrator of the form (|35 j) can be checked 
by comparing the series of differential operators ^(h) with the series exp(hF) 
associated to the flow of the system (fl]). That is, the integrator (|35[) is of order 
r if 

#n = —,F n , l<n<r. (44) 

Instead of using (|42|) to obtain the terms ^f n of the series *&(h), one can equiv- 
alently consider the formal equality 

Vj/^) _ e -Y(-hai) e Y(ha 2 ) . . . Q -Y(-ha 2a -i) ^Y{ha 2s ) ^ (45) 

to obtain the series expansion of log(^(/i)) = ^ n>1 h n F n , so that rth order 
compositions methods can also be characterized by the conditions 

Fi = F, F n = for 2 < n < r. (46) 

As for the splitting integrator ([36]) when the ODE ([I]) is split in two parts, 

f(x) = fW(x) + fW(x), (47) 

let and F^ be the Lie derivatives corresponding to and respectively, 
that is, 

p [a] [9] (*) = E /] al (*) ^ (*) • [5] (*) = E /j 61 (*) ^ (*) ( 48 ) 

for each 51 G C^M^IR) and each x G M D . Then, the series ^(h) of differential 
operators associated to the integrator iph in ()36[) can be formally written as 



3 Order conditions of splitting and composition meth- 
ods 

There are several procedures to get the order conditions for the coefficients of 
splitting and composition methods of a given order. These are, generally speak- 
ing, large systems of polynomial equations in the coefficients which are obtained 
from equations ([46 p . Perhaps the two most popular are (i) the expansion of 
the series log(^(h)) = Yln>ih n F n of vector fields by applying recursively the 
Baker-Campbell-Hausdorff (BCH) formula [88], and (ii) a generalization of the 
theory of rooted trees used in the theory of Runge-Kutta methods, which al- 
lows one to get an equivalent set of simpler order conditions in a systematic 
way [63] (see also [S])- In this section we first summarize briefly how to get 
these equations with the BCH formula, and then we present a novel approach, 
related to that in [63], but based on Lyndon words instead of rooted trees. 
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3.1 Order conditions via BCH formula 

As is well known, if X and Y are two non-commuting operators, the BCH 
formula establishes that formally, e x e Y = e z , where Z belongs to the Lie 
algebra C(X, Y) generated by X and Y with the commutator [X, Y] = XY — 
YX as Lie bracket [84] . Moreover, 

oo 

Z = log(e x e y ) =X + Y+ Z m , (50) 

m=2 

with Z m (X, Y) a homogeneous Lie polynomial in X and Y of degree m, i.e., it is 
a Q-linear combination of commutators of the form [Vi, [V2, • • • , [V m _i, V m ] . . .]] 
with Vi £ {X, Y} for 1 < i < m. The first terms read 

Z2 = \[X,Y] 

Z3 = ~[[X,Y],X\ + ±[[X,Y],Y] 
Za = ±[[[X,Y},Y],X], 

and explicit expressions up to m = 20 have been recently computed in an 
arbitrary generalized Hall basis of C(X, Y) [22]. 

The procedure to get the order conditions for the composition method (|35p 
with this approach can be summarized as follows. First, consider the series 
of differential operators ^>{h) associated to the integrator (j35|) . expressed as a 
product of exponentials of vector fields, i.e., equation ([43]) . Then, apply repeat- 
edly the BCH formula to get the series expansion \og{^{h)) = Yl n >l h n F n . In 
this way, one gets 

log(¥(/i)) = hw 1 Y 1 + h 2 W2Y 2 + h 3 (w 3 Y 3 + w 1 2[Y 1 ,Y 2 ]) (51) 

+h 4 ( W4 Y 4 + wi3\¥i,Y 3 ] + loiia [n, [ri,y 2 ]]) + 0(h 5 ) 

where the Wj v ..j m are polynomials of degree n = ji + ■ ■ ■ +j m in the parameters 
ai,... , a 2s- The first such polynomials are 

2s 2s 2s 

wi = J2 ai > ^2 = ^(-l) i a-, w 3 = J2 a l ( 52 ) 

i=l i=l i=l 

In general, the expressions of the polynomials Wj 1 --j m in (|51|) obtained by re- 
peated application of the BCH formula are rather cumbersome. 

The order conditions for the composition integrator (|35|) are then obtained 
by imposing equations (|46p to guarantee that the scheme has order r > 1. Thus, 
the order conditions are w\ = 1, and Wj 1 ...j m = whenever 2 < + - ■ ■+j m < r - 

One can proceed similarly to get the order conditions of the splitting scheme 
(|36p in terms of the coefficients aj, b^. Consider the series of differential operators 
^f(h) associated to the integrator (|36|) expressed as ([4^]) ; then, apply repeatedly 
the BCH formula to get the series expansion log(^(/i)) = X^n>i h n F n , so that 



15 



the order conditions will be obtained by imposing equations (|46p to guarantee 
order r > 1. It can be seen that the following holds for \og{%>{h)), 

log(tt(fc)) = ^(^FW + w^M) + li 2 ^FH + ^(t^F^ + v aba F^) 

+h\v abbb F^ + v abba F^ + v abaa F^) + 0(h B ), (53) 

where 

p[ab] _ [pM^Ml ^[066] _ [p[a6]^[&]i ^[a6o] _ [p[a&]^[a]i 
^[aftfeb] _ [p[a&&] ^[afeba] _ ^p[abb] ^ p[ah ^ p[abaa] _ rp[aba] ^ p[a]] 

and v a ,v b , v ab , v abb , v aba , v abbb , . . . are polynomials in the parameters a^, 6j of the 
splitting scheme ([36]) . In particular, one gets 

s s+1 



yZ a i> V b = '^2b i , V ab = -V a V b - ^ b i a j' ( 54 ) 

i=l i=l !<*<i<s 

~\ V a V b+ ^2 a i b j a ki 2v abb = -v a v b - ^ hajb k . 

l<i<j<k<s l<i<j<k<s+l 



2v aba 



From ()53p , we see that a characterization of the order of the splitting scheme 
(l36l) can be obtained by considering v a = v b = 1 and t> a & = v abb = v aba = • • • = 
up to polynomials of that form of the required order. The set of order con- 
ditions thus obtained will be independent in the general case if the vector 
fields F^\F^\F^ ab \F^ bb \F^ considered in ([53]) correspond to a basis of 
the free Lie algebra on the alphabet {a,b}. Notice that in (|53p we have con- 
sidered a Hall basis (the classical basis of P. Hall) associated to the Hall words 
a, b, ab, abb, aba, abbb, abba, abaa, ■ ■ ■ [69]. The coefficients v w in (|53p correspond- 
ing to each Hall word w can be systematically obtained using the results in [62] 
in terms of rooted trees and iterated integrals. An efficient algorithm (based 
on the results in |62j) of the BCH formula and related calculations that allows 
one to obtain (|53p up to terms of arbitrarily high degree is presented in [22]. 



3.2 A set of independent order conditions 

We next present a set of order conditions for composition integrators (135 1) de- 
rived in [25j . 

From ([4l])-([43]), it follows that 

V(h) = I + J2h n E u jl ... jr (a 1 ,...,a 2s )X jl ---X jr , (55) 

ri>l jiH Yjr=n 

for some polynomial functions Uj x ...j r of the parameters a\, . . . , a 2s of the method 
We next introduce some notation in order to explicitly write these polynomials. 
For each positive integer j, we write j* = j — 1 if j is even, and j* = j if j is odd. 
Finally, for each pair (i, j) of positive integers, we write = (— l)^ % ~ l \aj) 1 . 

That is, (Xj = {o.j) % if j is even or i is odd, and = — {a.j) % if j is odd and i 
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is even. Now, it is not difficult to check that, for each multi-index (ii, . . . , i m ) 
of length m > 1 and (a±, . . . , a 2s ) £ M. 2s , 

u il ... im (a 1) ...,a 2s )= J2 4i ) -'- a t ) - ( 56 ) 

1 <jl <3{ <h < ■ ■ ■ < 3m - 1 <3* m _ X <jm < 2s 

Obviously, each u^...^ can be seen as a real-valued function defined on the set 

{(ai,...,a 2a ) £l 2s : s > 1}. (57) 

Observe that each Uj r ..j m (ai, . . . ,a 2s ) is a polynomial of degree ii — ii -\- • • • ~\~irn 
in the variables ai, . . . , «2s- 

Now, the order conditions of the composition scheme (|35[) can be obtained 
by comparing the series (|55j) with exp(/ii ? ), that is, (|44l) . Since Xi = i 7 , as the 
basic integrator x/i is assumed to be of order 1, we have that the method is of 
order r if for each multi-index . . . , i m ) with ii + • • • + i m = n < r, 

u il ..^ m (ai J ...,a 2a )-| Q otherwige . (58) 

However, such order conditions are not independent. For instance, it can be 
checked that 

l r 2 x 1 3 1 1 

U H = o( u l + u 2), U21 = -U12 + U 3 + U±U 2 , UlU = -«i + -U12 + ~U S , 

2 bid 

which implies that the order conditions ([35]) for tin, U12, iim are fulfilled 
provided that the conditions for u\,va,uz, u±2 hold. 

A set of independent order conditions can be obtained as follows. Consider 
the lexicographical order < (i.e., the order used when ordering words in the 
dictionary) on the set of multi- indices. A multi-index . . . ,i m ) is a Lyndon 
multi-index if . . . , i^) < (ik+i, • • • , i m ) f° r each 1 < k < m. For each n > 1, 
we denote as L n the set of functions Uj r ..j m such that . . . ,i m ) is a Lyndon 
multi-index satisfying that i\ + ■ • • + i m = n. The first sets L n are 

L\ = {u±}, L 2 = {u 2 }, L 3 = {u 12 ,u 3 }, L 4 = {u U2 ,u 13 ,U4}, 
L5 = {uui 2 ,un 3 ,ui 22 ,u 14i ,u 23 ,u 5 }. 
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In particular, we have 

ui(a>i, ... , a 2s ) 

u 2 (a>i, . . . ,a 2s ) 

u 3 (ai, . . . ,a 2s ) 

ui 2 (a>i, . . . ,a 2s ) 

u^ai, . . . ,a 2s ) 



u l3 (ai, . . . ,a 2s ) 

nn 2 (ai, . . . ,a 2s ) 
and so on. 

We can finally state the following result [25J : Given (ai, . . . , a 2s ), the inte- 
grator (|35p is of order r for arbitrary ODE systems (JTJ) and arbitrary consistent 
integrators \h if and only if a,\ + • • • + a 2s = 1 (i.e. ui(a±, . . . , a 2s ) = 1) and 

r 

Vii£|jL n , u(ai, . . . ,a 2s ) = 0. (59) 

n>2 

Furthermore, such order conditions are independent to each other if arbitrary 
sequences (a±, . . . , a 2s ) of coefficients of the method are considered. 

3.3 Order conditions of compositions methods with symmetry 

The order conditions are simplified for (2s)-tuplas (a±, . . . ,a 2s ) such that 

a 2s -j+i = aj, for all j. (60) 



2s 




3=1 




2s 








3=1 




2s 




rj 

/ / a 3 > 




3=1 




2s 32 




/ , ^ ^ "32 ^ 


fX ■ 

u 3i ' 


32=1 j!=l 












3=1 




2s 3'| 




E a 3 3 2 E°^' 




32=1 ji = l 




2s 33 


32* 






33=1 j 2 = l 


31 



It is easy to check that the simplifying assumption (|60j) implies that the com- 
position integrator (|35p is time-symmetric (i.e., ip^ = iph)- In that case, only 
the conditions for u £ L n with odd n remain independent. 

The order conditions can be alternatively simplified by requiring that 

a 2 j = a 2j -i, Vj, (61) 

in which case, only the conditions for Lyndon multi- indices (ix, . . . ,i m ) with 
odd ii, . . . ,i m are required. The simplifying assumption (|6ip means that the 
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k 


1 


2 


3 


4 


5 


6 


7 


8 


9 


10 


11 




1 


1 


2 


3 


6 


9 


18 


30 


56 


99 


186 


m k 


1 





1 


1 


2 


2 


4 


5 


8 


11 


17 



Table 1: The numbers n k and m k of independent order conditions for general 
composition methods (|35p and for compositions (|62[) of a basic time-symmetric 
method, respectively. 



composition integrator (|35p can be rewritten as 

where /3j = 2«2j and 5^ is the self-adjoint second order integrator S h = 

Xh/2°X* h/2 - 

The order conditions are thus considerably reduced if one considers compo- 
sition methods satisfying both assumptions (|60p - (|6ip . that is, methods of the 
form (|62p satisfying that 

Pj = p a - j+1 , Vj. (63) 

Schemes of this form can be dubbed as symmetric compositions of symmetric 
schemes. For instance, for order r > 6 one has the conditions 

2s 2s 2s 2s 33 32 

Y a 3 = Y ° 3 3 = °> Y a ) = °' Y a % Y a n Y a n = 

3=1 3=1 3=1 33=1 ,j2=l 31 

in terms of the a% coefficients (the actual expressions in terms of (3% are slightly 
more involved). In Table Q] we display for each k > 1 the number n k of Lyndon 
multi-indices (ii, . . . , i m ) with «i + • • • + i m = k, and the number of Lyndon 
multi-indices . . . , i m ) with + - ■ ■+i m = k and odd indices i\, . . . , i m . Thus, 
the number of independent conditions to guarantee that the general composition 
integrator (|35p is at least of order r is n\ + • • • + n r , while in the case of the com- 
position (162p based on a symmetric second order integrator S h (or equivalently, 
a composition integrator (|35p with the additional symmetry condition (|6ip ). the 
number of independent order conditions is mi + • • • + m r . If time-symmetry 
is imposed in the method (|35p (resp. (|62p ) by the additional assumption (|60p 
(resp. ([63]) ). then there are n\ + 713 + • • • + n2/_i (resp. m\ + 7713 + • • • + m,2Z-i) 
independent conditions that guarantee order at least r = 21. 



3.4 Relation among different sets of order conditions of com- 
position methods 

In [63], a set of independent necessary and sufficient order conditions is given 
using labelled rooted trees (see also [41]). A family of sets {'7^}n=i,2,... of func- 
tions defined on the set (|57|) is identified such that the integrator (|35|) is of order 
r if and only if a± + ■ ■ ■ + Q2 S = 1 and 

r 

Vue(j7^, u(a t , . . . ,a 2s ) = 0. (64) 

n>2 
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Each u(a±, . . . , a 2s ) with u G T n is (as in the case where u G £ n ) 5 a polynomial 
of homogeneous degree n. In particular, 

Ti = {-ui}, T 2 = {-u 2 }, T 3 = {ii 2 i,n 3 }, T A = {t) 2 n, u 3 i,u 4 }, 

where the functions of the form Ui t ...i m are defined in ([56]) . and 



As shown in [25] , the order conditions are equivalent to the conditions ([59]) , 
as both U n >iL n and U^i?^ generate the same graded algebra Ti = © n>1 Ti n 
of functions on the set ([57]) (for each u G 7i n , it(ai, . . . , «2s) is a polynomial of 
homogeneous degree n, actually, a linear combination of polynomials u^...^ of 
homegeneous degre n = %\ + • • • + i m ). For instance, it can be seen that 

v 2U = 2u 2U - u 22 = 2(n 112 - ui 3 + uiu X2 + u 3 ui) + u\u 2 + ^(u 4 - u^). 

Finding an independent set of order conditions for composition integrators is 
equivalent to finding a set of functions of homogeneous degree that generate 
the algebra Ti (see |25j ) for more details. 

Of course, the functions w^...^ in (|5ip obtained when deriving the order 
conditions of composition integrators by repeated use of the BCH formula also 
belong to the same algebra of functions. For instance, w n = u n , and wi 2 = 

Ul2 -U 3 - U X U 2 . 

Recall that Theorem [T] characterizes the order conditions of splitting in- 
tegrators of the form (|36p . where the ODE ([1]) is split in two parts (|47p . in 
terms of the order conditions of composition integrators (I35p . Actually, the 
polynomials v a ,v b , v ba ,v baa , v bab , v baaa , ... (on the parameters a^bi) in ([53]) can 
be rewritten as linear combinations of the polynomials (on the parameters ctj) 
in (]56p provided that (|37p and v a = v b = u\ hold. In particular, it can be seen 
that 

u 2 

v a b = y, 

Vabb = JX (-U 3 - 3ni2 + 3«2l) , 

1 . 

Vaba = — {U 3 ~ 3u U + 3u 21 ) , 

1 , . 

Vabbb = — \U 22 - U 3 i + U U2 - 2U121 + U 2n ) , 

Vabba = 24 (~ Ui ~ 2Ul3 + 4U22 ~ 2u ' il + 4ni12 ~ 8Ul21 + 4u 21l) , 
1 , 

Vabaa = {-U13 + U 22 + «112 ~ 2^121 + «21l) ■ 

3.5 Negative time steps 

It has been noticed that some of the coefficients in splitting schemes (|36p are 
negative when the order r > 3. In other words, the methods always involve 
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stepping backwards in time. This constitutes a problem when the differential 
equation is defined in a semigroup, as arises sometimes in applications, since 
then the method can only be conditionally stable [57]. Also schemes with 
negative coefficients may not be well-posed when applied to PDEs involving 
unbounded operators. 

The existence of backward fractional time steps in this class of methods 
is unavoidable, as shown in [35j [75j [79]. In fact, it can be established in an 
elementary way by virtue of the relationship between the order conditions of 
schemes f)36f) and (I35p stated in Theorem [T] [4j: Any splitting method of the 
form (|36p that has order r > 3 neccesarily must fullfil the condition 

2s s 

u 3 (ai, . . . , a 2s ) = ^2 of = ^(ali-i + o? 2i ) = 0, (65) 

8=1 1=1 

with coefficients ctj obtained from the relations (|37p . Since, for all x, y £ K, it is 
true that x 3 + y 3 < implies x+y < 0, then there must exist some i 6 {1, . . . , s} 
in the sum of (J65j such that 

a 2i _ l + a 2i < and thus o-n-i + «2i = a* < 0. 

Obviously, one can also write (by taking oq = 0) 

2s+l s+1 

u 3 (ai, . . . , a 2s ) = ^ a? = 5^(oii-i + "2^-2) = 
i=0 i=l 

just by grouping terms in a different way, and thus, by repeating the argument, 
there must exist some j € {1, . . . , s + 1} such that 

a 2 j-i + a 2 j-2 = bj < 0. 

This proof shows clearly the origin of the existence of backward time steps: the 
equation U3 = can be satisfied only if at least one and one b{ are negative. 
According to this conclusion, any splitting method of the form (|36p verifying 
the order condition 113 = has necessarily some negative coefficient and also 
some negative b{. 

3.6 Near-integrable systems 

In Hamiltonian dynamics one often encounters systems whose Hamiltonian 
function H is a small perturbation of an exactly integrable Hamiltonian Hq, 
that is H = Hq + eH\ with e« 1. The perturbed Kepler problem (i20|) belongs 
to this category of near-integrable Hamiltonian systems. The gravitational N- 
body problem (|2ip . when using Jacobi coordinates, also falls within this class 
of problems. In that case, Hq represents the Keplerian motion and eH\ the 
mutual perturbations of the bodies on one another [86]. 
More generally, let us consider an ODE system 

s' = /M(s) + e/M(x), (66) 
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containing a small parameter |e| <C 1. If the exact /i-flows cpl ' and <p\' of 
x > _ y[ a ]( x ) anc l x ' = e/[ fe l(x) respectively can be efficiently computed, then 
a scheme tph of the form (|36|) can perform particularly well provided that the 
coefficients ai,bi are appropriately chosen. To see this, consider the Lie deriva- 
tives (|48p of 7^ and respectively, so that the corresponding series ^(/t) 
(|49p of differential operators associated to the scheme (|36p becomes 

= e biheFW e aihFW , _ , e b s heF^ & a s hFM e b s+1 heF^ 

Successive application of the BCH formula then leads to (|53p with replaced 
by eF®, that is 

log(*(/ i )) = + e(^F [6] + tfiW^ + ft 3 ^ a fW + h 4 v abaa F^) 

In practical applications one usually has e <C /i (or at least e ~ /i), so that 
one is mainly interested in eliminating error terms with small powers of e. For 
instance, if the coefficients a^&j of the splitting methods are chosen in such a 
way that 

v a = l = v b , v ab = v aba = v abaa = v abb = 0, 

then 

log(*(») -hF = 0{eh 5 + e 2 h 4 ), 

where F = F^ + eF*-". More generally, one is interested in designing methods 
such that [12] 

log(*(h)) -hF = 0{eh Sl+1 + e 2 h S2+1 + e 3 h S3+1 + ■■■ + e m h s ™ +l ). (67) 

Observe that s\ is the order of consistency the method would have in the limit 
e — > 0. It is relatively easy to eliminate errors of order eh k because there is 
only one such term for each order k, namely h k e v aba ... a F^ aba "' a ^ (with Fl aba — a \ = 
[[...[[FW,F%fW]...],fW]). 

If one is interested in designing methods that approximate the exact solution 
up to higher powers of e, more terms have to be considered. In particular, there 
are |_^(fc — 1)J terms of order 0(e 2 h k ) and [g(A; — 1)(A; — 2)J terms of order 
0(e 3 h k ) [53]. 

3.7 Runge— Kutta— Nystrom methods 

Suppose now that one is interested in integrating numerically second-order ODE 
systems of the form 

y" = g(y), (68) 

where y G R d and g : R d — > R d . In this case it is still possible to use schemes 
(|36p applied to the equivalent first-order ODE system. More specifically, intro- 
ducing the new variables x = (y, v), with v' = y, and the maps : M. 2d — > M. 2d 
and /M : R 2d -» M. 2d defined as 

fW(v,v) = (v,0), fW(y,v) = (0,g(y)), (69) 
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k 


2 


3 


4 


5 


6 


7 


8 


9 


10 


11 




1 


2 


3 


6 


9 


18 


30 


56 


99 


186 


h 


1 


2 


2 


4 


5 


10 


14 


25 


39 


69 



Table 2: The numbers Ik of independent order conditions (at order k) for split- 
ting methods in the RKN case compared to the numbers n\~ of conditions in 
the general case. 

equation (1681) can be rewritten 

' = f [a] {x) + f [b] (x). Then, the splitting 
scheme (|36p can be efficiently implemented, as the exact h-Rows cp^ and ipf? 
of x' = f^(x) and x' = p b *(x) are simply given by 

<P l h(y,v) = (y + hv,v), 

<p [ h(y> v ) = (Vi v + h 9{y))- 

It is not difficult to check that the splitting schemes of the form (1361) 
are particular instances of Runge-Kutta-Nystrdm (RKN) methods (see for in- 
stance [S]). 

One of the most important applications of this class of schemes is the study 
of Hamiltonian systems of the form H(q,p) = T(p) + V(q), where the kinetic 
energy T(jp) is quadratic in the momenta p, i.e., T(p) = ^p T Mp for a sym- 
metric square constant matrix M, and V(q) is the potential. In that case, the 
corresponding Hamiltonian system can be written in the form f|68j) with y = q, 
y' = v = Mp, and g(y) = -W(y). 

Although a splitting integrator (I36p designed for arbitrary ODE systems 
x' = f(x) split into two parts (|47p will perform well when applied to a second 
order ODE system of the form (I68p with the splitting (169)) . much more efficient 
methods can be designed in that case [56)119]. The main point here is that in the 
case UBS) (we will refer to that as the RKN case), [[[F® , F^}, F®], F®] = 
identically. This is equivalent to F^ habb ^ = in (|53p . which introduces some 
linear dependencies among higher order terms in the expansion of \og(^(h)) 
(see |59] for a detailed study). This means that the characterization given in 
Theorem [T] for a splitting integrator ([36]) to be of order r (for r > 4) is no longer 
applicable if one restricts to the case (|69|) . In Tabled! the number of necessary 
and sufficient independent order conditions for a splitting method (|36|) to be of 
order r in the RKN case is compared to the general case: For arbitrary systems 
split into two parts (|4T|) . there are 2 + ri2 + • • • + n r independent conditions 
(including the two consistency conditions v a = u& = 1), while in the RKN case, 
the number of independent order conditions is 2 + li + • • • + l r . Unfortunately, 
up to order three the order conditions are the same in both cases and then, the 
results for negative time steps still apply. 

Since the reduction in the number of order conditions is due to the fact 
that f^ a \y, v) is linear in v, it is immediate to see that the methods obtained 
in this way also apply to the more general problem f^ a \y,v) = (wi(x)v + 
W2(y), wz{x)v + Wiiyj), which includes the system 

y" = My' + g 1 (y)+g 2 (y). (71) 
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Here splitting RKN methods are useful if the reduced problem y" = My' +g\(y) 
(i.e., f^ a \y,v) = (v,Mv + gi(y))) is easily solvable. For Hamiltonian systems, 
this generalization corresponds to H(q,p) = T(q,p) + V(q), where T(q,p) = 
^p T M(q)p + f T (q)p + W(q). Obviously, the exact solution for T(q,p) is only 
known for some particular cases, e.g. if T corresponds to the Kepler problem 
in (f20|) or (|2Tj) or to the harmonic oscillator in (fT9|) or (f22j) . 

It is interesting to note that in quantum mechanics the kinetic and potential 
energy verify analogue commutator rules to classical mechanics and then RKN 
methods can also be used. If applied, for instance to problems (f26|) and (p7|) . 
one should keep in mind that in the resulting composition method ipff must 
correspond to the kinetic part. 

We should also remark that, if the Hamiltonian function H(p, q) = ^p T Mp+ 
V(q) is such that in addition V(q) = ^q T Nq (i.e., it corresponds to the gener- 
alized harmonic oscillator (|18|) ). then the number of order conditions reduces 
drastically: it is not difficult to see that there is only one independent condition 
to increase the order from r = 2k — 1 to r = 2k, and only two to increase the 
order from r = 2k to r = 2k + 1 (see [9] for more details). We will return later 
to this system and take profit of its special features to design specially adapted 
splitting methods. 

4 Additional techniques to reduce the number of or- 
der conditions 

4.1 Methods with modified potentials 

The splitting method ([36]) can be generalized by composing the exact flows 
of other vector fields in addition to and provided that they lie on 
the Lie algebra generated by and F^. For instance, one could consider 
compositions that, in addition to </?j^ and ip^\ use the /i-flow of the vector 
field F^ = [[FW,fI%FI% To illustrate this fact, consider the composition 

/ [b] [a] [b] [abb] [b] [a] [b] f „ n \ 

The scheme ([72]) . constructed in [261 H7], is °f order four. Indeed, by repeated 
application of the BCH formula to 

T/ ,x hplt] h F [ a ] h F [b] _jLp[abb] hp[b] hp[a] h p[b] 

ty(h) = e$ ea es e ™ es e2 e« 

one can check that *(h) = e h ^+^)+0(h^ _ 

Recall that, although the h- flows cp^ and tp^ of the vector fields F^ and 
FW are by assumption computed easily, this is not necessarily the case for the 
h-Row </?j^ b6 ' of the vector field F^ abb ^ . However, in the RKN case (|68p considered 
in Subsection E21 where / [a] and f [b] are of the form JES}, the h-ftow of F^ 
is of the form 

^ bb \y, v) = (y,v + h 3 gW(y)), where g®(y) = 2g'{y)g{y). 
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This shows in addition that and (^Sf 6 ^ commute, and that in particular, for 
arbitrary bj,cj £ R, 

$h ° ° . «) = (f . « + 2/lb i + ^(VMV)), (73) 



which is precisely the /i-flow of the vector field 2b j + Cjh 2 F^ abb \ It thus 
makes sense to construct methods defined as compositions of \ and maps of 
the form ([75]) for j = 1, . . . , s. 

For Hamiltonian systems (?) = r(p) + V(q) with quadratic kinetic 
energy T(p) = \p T Mp, the vector field = [[F^,F^],F^] is the vector 

field associated to the Hamiltonian function {SJV) T W, which only depends 
on the position vector q. Thus, (|73p is just the /i-flow of the system with 
Hamiltonian function 

2b, V(q) + Cj h 2 (VV(q)) T VV(q), (74) 

which reduces to the potential V(q) of the system for bj = 1/2 and Cj = 0. 
This explains the term 'splitting methods with modified potentials' used in the 
recent literature [5T] [7_U [87] to refer to splitting methods obtained by composing 
the /i-flows of T and modified potentials of the form (|74|) . 

This procedure can be generalized by considering "modified potentials" of 
higher degree in h. In particular, the flow l p^ bbab ^ Q f the vector field 

p[abbab] _ mpW } jrV>]] p[bh \F^ a \F^]], (75) 
is of the form (p^ bbab ^ (y 5 v) = (y,v + h 5 g^ (y)), and similarly for the vector fields 

p[abbabab] _ [[[[_p[o] ? pV>] j ^ p[b]^ ^[a] ^ ^[6]jj J^N ^ _p[6]]j ^ (75) 

with /i 5 g[ 5 '(y) replaced by h 7 g^ 7 ^(y) and /i 7 ^ 7 ' 1 ! (y) respectively. The functions 
be written in terms of g and its partial derivatives (see [13] 

for more details). 

In some applications, the simultaneous evaluation of g(y), g^(y), g^(y), 
9^ 7 ' l \y) an d g^ ,2 \y) is not substantially more expensive in terms of computa- 
tional cost than the evaluation of g(y) alone. In that case, by replacing in the 
scheme ([36]) each tpf^ h by the /i-flow of 



C h (bi, a, di, en, ea) = h F® + h 2 Cl F^ + h% 

+h 6 (e iA F [abbabab] + e i>2 F [abbaabb] ) (77) 

additional free parameters are introduced to the scheme without increasing too 
much the computational cost, which allows the construction of more efficient 
integrators. 

Of course, this can be further generalized by considering more general nested 
commmutators of F^ and F^ that gives rise to "modified potentials". In that 
case, higher degree commutators afected by higher powers of h should be added 
in JZ7J). 
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Notice that in this case the coefficients aj, 6, have not to satisfy all the order 
conditions at order r > 3 and then, the results for negative time steps do not 
apply in this case. As a result, schemes with positive coefficients do exist. In 
this case, negative coefficients appear in methods of order six [27] , 



4.2 Methods with processing 

Recently, the processing technique has been used to find composition meth- 
ods requiring less evaluations than conventional schemes of order r. The idea 
consists in enhancing an integrator iph (the kernel) with a parametric map 
7T/j : MP — ► R D (the post-processor) as 

= n h o ?p h o vr" 1 . (78) 

Application of n steps of the new (and hopefully better) integrator ip^ leads to 

which can be considered as a /i-dependent change of coordinates in phase space. 
Observe that processing is advantageous if iph is a more accurate method than 
ipfi and, either the cost of 7r/j is negligible or frequent output is not required, 
since in that case, it provides the accuracy of iph at essentially the cost of the 
least accurate method 

The simplest example of a processed integrator is provided in fact by the 
Stormer-Verlet method. As a consequence of the group property of the exact 
flow, we have 

c [2] [a] [b] [a] [a] [b] [a] [a] [a] 

S h = <Ph/2 0( Ph 0( Ph/2 = ( Ph/2 0( Ph 0( Ph °<P-h 0< Ph/2 

[a] [a] -1 fr7n \ 

= V>h/2 <p_ h/2 =K h °Xhon h (79) 

with 7T/j = <^H 2 and the symplectic Euler method \h = <P$- Hence, 

applying the basic integrator Xh = <P$ ° P$ with processing yields a second 
order of approximation. 

Although initially proposed for Runge-Kutta methods p2], the processing 
technique has proved its usefulness mainly in the context of geometric numerical 
integration [41j . where constant step-sizes are widely employed. 

We say that the method iph is of effective order r if a post-processor 717,, 
exists for which iph is of (conventional) order r [18] . that is, 

tt^ o o TT" 1 = (p h + 0(h r+1 ). 

Hence, as the previous example shows, the basic splitting <p$°<p$ is of effective 
order 2. Obviously, a method of order r is also of effective order r (taking 
7T/j = id) or higher, but the converse is not true in general. 

The analysis of the order conditions of the method iph shows that many 
of them can be satisfied by ir^, so that iph must fulfill a much reduced set of 
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restrictions [61 [TT]. For instance, if the kernel is defined as (|35p with a basic 
first order integrator Xh and the post-processor is similarly defined as 



*h = X l2m h o X* 2s _ lh ■ ■ ■ X J2 h ° X* lh 



(80) 



then, conditions 



ii\{a) = 1, 112(a) = 113(a) = 114(a) = 



(81) 



guarantee that the kernel iph is of effective order four. If in addition the post- 
processor ([80|) satisfies 

111(7) = 0, 112(7) = "12(a), 1*3(7) = "13(a), "12(7) = ""112(a), 

then the processed integrator (|78p has conventional order four. Here, we use 
the notation a = (a\, . . . ,a2 S ) and 7 = (71,..., 72m) for the coefficients of 
the kernel and the post-processor respectively. If in addition the following 
conditions are fulfilled by the coefficients of the kernel, 



then the kernel iph has at least effective order five. In that case, the processed 
method (178p achieves conventional order five if in addition, the equalities 



hold for the coefficients of the post-processor (|80l) . 

Thus, the number and complexity of the conditions to be verified by the 
coefficients ay of a kernel of the form (I35j) is notably reduced. Highly efficient 
processed composition methods that take advantage of that have been proposed 
[111 [57] . Nevertheless, when both the kernel iph and the post-processor ir^ are 
constructed as a composition of the form ([35]) (or ([36]) ). the use of the resulting 
processed scheme is not recommended in situations where intermediate results 
are required at each step. Indeed, the total number of compositions per step in a 
processed method (178p of that form is typically higher than for a non-processed 
method of comparable accuracy. 

To overcome this drawback, in [6] a technique has been developed for ob- 
taining approximations to the post-processor at virtually cost free and without 
loss of accuracy. The key idea is to replace tt^ by a new map tt^ ~ 7T/j obtained 
from the intermediate stages in the computation of iph- The post-processor 
can safely be replaced by an approximation tt^, since the error introduced by 
the cheap approximation tt^ is of a purely local nature [6] (it is not propagated 
along the evolution, contrarily to the error in tt^ 1 )- 

In [6], a general study of the number of independent effective order order 
conditions versus the number of conventional order conditions is presented. In 
particular, it is shown that, in the case of kernels of the form ()35p . the number 
of conditions to increase the effective order of the kernel from k > 1 (resp. 



u 5 (a) = M 2 3(a) = 0, 



2«i22(a) + uu(a) + uu{a) 2 = 0, 



114(7) = "14(a), 1113(7) = ""113(a), 
"112(7) = "1112(a) + -ui 2 (a) 2 - -«n 2 (a) 



^"12(0) 
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r 


2 


4 


6 


8 


10 


12 


N r 


1 


3 


9 


27 


83 


269 


N* 


1 


2 


5 


14 


40 


127 


M r 


1 


2 


4 


8 


16 


33 


M* 


1 


2 


3 


5 


8 


14 


L r 


2 


4 


8 


18 


43 


112 


L* 


2 


3 


5 


10 


21 


51 



Table 3: Number of conventional and effective order conditions for symmetric 
kernels: (i) (N r ,N*) for composition integrators (|35|). (ii) (M r ,M*) for compo- 
sitions (|62p of symmetric second order basic integrator, (iii) (L r , L*) for splitting 
integrators in the RKN case. 

k = 1) to k + 1 is rik + i — (resp. r&2 — rai + 1), where each is the cardinal 
of Lfc, that is, the number Lyndon multi- indices of degree k. Thus, whereas 
the total number of independent conditions to achieve conventional order r is 
n\ + • • • + n r , only 1 + n r conditions have to be imposed to the kernel for 
effective order r. If the kernel (|35j) is time-symmetric (i.e., if its coefficients 
satisfy ([60]) ). then there are N r = Yli=i n 2j-i independent conditions for order 
r = 2q, and N* = n\ + Ylt=i ( n 2j+i — n 2j) conditions for effective order r = 2q. 
A similar situation occurs for the total numbers M r and M* of conventional 
and effective order conditions of symmetric kernels of the form (|62p with (|63[) 
(where the are replaced by the number mt in Table [T]). That also happens 
to be true for symmetric kernels of the form (|36p . both in the general case 
(which is essentially equivalent to the case of kernels of the form (|35p ) and in 
the RKN case considered in Subsection 13.71 In Table El the total number of 
conditions for conventional order r = 2q for symmetric kernels is compared with 
the total number of effective order conditions in three kinds of integrators: (i) 
(N r , N*) for composition (|35p of a basic first order integrator and its adjoint, (ii) 
(M r ,M*) for compositions (|62p of a symmetric second order basic integrator, 
(iii) (L r ,L*) for splitting integrators in the RKN case. 

5 A collection of splitting methods 

As we have mentioned before, splitting methods have found application in many 
different areas of science during the last decades. It is therefore not surprising 
that there is a large number of different schemes available in the literature. 
Sometimes, even the same method has been rediscovered several times in dif- 
ferent contexts. Our aim in this section is to offer the reader a comprehensive 
overview of the existing methods, by classifying them into different families and 
giving the appropriate references where the corresponding coefficients can be 
found. 

At this point it is important to remark that the efficiency of a method is 
measured by taking into account the computational cost required to achieve 
a given accuracy (we do not take into account the important property of the 
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stability of the methods). For instance, given several methods of order r with 
different computational cost (usually measured as the number of stages or evalu- 
ations of the functions involved) , the most efficient method does not necessarily 
correspond to the cheapest method. The extra cost of some methods can be 
compensated by an improvement in the accuracy obtained. 

We next present a short review indicating the splitting methods which have 
been published in the literature at different orders, with different number of 
stages and for several families of problems. 

Symmetric compositions of symmetric methods. As we pointed out in 
section 12. 1| although by applying recursively the composition f|30|) - (|31|) it is 
possible to increase the order, the resulting methods are computationally ex- 
pensive. To reduce the number of evaluations the more general composition 
(|62p may be considered to achieve a given order r. If we choose symmetric 
compositions {(3 s +i-i = Pi), then half of the parameters of the method are 
fixed, but the order conditions at even orders are automatically satisfied. In 
other words, the parameters of a (non-symmetric) method of order r = 2k have 
to solve a system of X^=i m « equations (see Table [3]), whereas for a symmet- 
ric composition e r = 777,2 + 777.4 + • • • + rri2k order conditions are automatically 
satisfied if the order conditions at odd orders are fulfilled. In this way, only 
M r = mi + 7713 + • • • + m 2 k-i independent order conditions need to be imposed 
in the case of symmetric compositions. Due to this fact, the number of con- 
ditions to be solved (which is typically the bottleneck in the numerical search 
of methods) is reduced considerably when imposing symmetry. Furthermore, 
since m 2 i~i < rri2i then M r < e r and symmetric compositions, in addition 
to having more favourable geometric properties (due to the time-symmetric 
property), usually require smaller number of stages than their non-symmetric 
counterparts. Taking into account the number M r (resp. M*) of independent 
conditions to achieve conventional order r (resp. effective order r) from Table El 
it is possible to determine the minimum number s r = 2M r — 1 of stages of the 
integrator (resp. the minimum number k r = 2M* — 1 of stages for the kernel) 
required by a method of order r (resp. effective order r) 

In this way one has to solve a system of M r or M* nonlinear polynomial 
equations with the same number of unknowns The number of real solutions 
typically increase a good deal with r. In general, these equations have to be 
solved numerically and getting all solutions is a very challenging task, even for 
moderate values of r. Once a number of solutions for the parameters Pi have 
been obtained, there remains to select that solution one expects will give the 
best performance when applied on practical problems, typically by minimizing 
some objective function. What is the most appropriate objective function in 
this case? A frequently used criterion is to choose the solution which minimizes 

If one takes additional stages in (|62p . for instance s = s r + 2, then one has 
an extra free parameter (notice the scheme is symmetric and two stages are 
required to introduce one parameter). By choosing Pi as this free parameter, 
then it is clear that 1-parameter families of solutions are obtained. For instance, 
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Order 


4 


6 


8 


10 


3-P3H2HISHIIS! 


7- [88] 


15-[H1 EU EH SS] 


31-dU [12 SUES] 




9- [MISS] 


17-IM1I35] 


33-03 SQlEaiZS] 




ll-13-[76] 


19-21-[76] 


35-001 [76] 






24-[20j 




P:3-17-[55j 


P:5-15-[7J 


P:9-19-I7] 


P:15-25-[7J 



Table 4: Symmetric compositions of symmetric methods published in the lit- 
erature. We indicate the number of stages (in boldface) and the pertinent 
reference. Processed methods are preceded by P. 



taking (3\ = one has the previous solutions and by continuation it is possible 
to get several of these 1-parameter families of solutions, but this procedure does 
not guarantee to find all solutions. 

Finally, one has to select that solution minimizing the value of C. Of course, 
additional stages can be introduced and the process is similar but technically 
much more involved. This objective function allows one to find very efficient 
methods involving additional stages, although the efficiency of methods with 
the same order but different number of stages cannot be compared from the 
value obtained for C. 

When we stop including additional stages in the composition ([62]) ? Two cri- 
teria are possible: (i) when one has enough stages available to achieve a higher 
order; (ii) when the performance of the actual methods constructed with addi- 
tional stages do not show a significant improvement in numerical experiments. 

For instance, the simple 4th-order scheme (|29h can be improved just by the 
5-stage generalized composition [78] 

c[2] „ o[2] q [2] q [2] q [2] , . 

with a = 1/(4 — 4 1//3 ), (3 = 1 — 4a, as numerical experiments clearly indicate. 
This is a particular case of (|62p where the value of C reaches a minimum and 
if we add two new stages with a new parameter then a 6th-order method can 
be obtained. 

In Table 0] we collect some of the most relevant methods from the literature 
with different orders and number of stages. At each order, r, we label the 
methods by the number of stages s and the reference where this method can be 
found. We also include methods obtained by using the processing technique, 
which are referred as P:s, where s is the number of stages for the kernel. We 
write S1-S2 for indicating that methods from si up to S2 stages are analyzed in 
that particular reference. 

Splitting into two parts. Composition of method and its adjoint. 

Next we review methods of the form (I36p (for ODEs that can be split into two 
parts) and (|35p . It is important to emphasize that, although the order condi- 
tions for both classes of methods are equivalent, the optimization procedures 
carried out to identify the most efficient schemes may differ. In consequence, 
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a particular method optimized for equations separable into two parts is not 
necessarily the best choice for a composition (|35p . although their performances 
are closely related. 

Considering, as before, symmetric compositions, i.e., a s+ \-i = a^, b s+ 2~i = 
bi in ([55]) and a 2s +i-i = on in (|3l)j) . it is easy to verify, from Table El that the 
minimum number of stages required to get a method of order r is s r = N r and 
of effective order r it is k r = N*. 

Note that schemes of order six or higher require more stages than compo- 
sitions (|62p . and only fourth-order methods seem promising. Nevertheless, one 
should recall that by including additional stages more efficient methods could 
be obtained. For instance, sixth-order methods require at least 9 stages (unless 
they are considered as composition of symmetric-symmetric methods in which 
case the 9 equations can be solved with only 7 unknowns) and the coefficients 
cti , bi or «j have to solve a system of eight nonlinear equations (in addition to 
consistency conditions). These equations have a very large number of solutions 
and it might be the case that one of them could correspond to a method with 
very small error terms. 

One optimization criterion frequently used when dealing with composition 
(f36|) is to work with the homogeneous subspace £ r +i = {F r +\ t \, . . . , F r +i jTlr+1 ) 
(where by -F r +i,i we denote the elements of the basis of the Lie algebra gen- 
erated by FW, F^ at order r + 1) and the leading error term, which can be 
expressed as YH=l c iF r +i,i- In this setting, one selects the solution minimizing 

E r+ \ = (X^I^i 1 l *! 2 ) 1 ^ 2 - This optimization criterion allows one to compare 

the performance of methods with different number of stages by introducing the 

1 It 

effective error, £f = sE r+1 , which normalizes with respect to the number of 
stages. 

For the composition <\35h . it is not so clear how to assign a weight to each 
element of the associated Lie algebra since their contribution on the error can 
differ significantly. One accepted choice consists in minimizing the objective 
function C = Y^aLi \ a i\- 

Methods up to order six built by applying this procedure can be found in the 
literature. They show for most problems a better efficiency than compositions 
(|62[) at the same order when applied to the same class of problems. We collect 
some of the most relevant schemes in Table [5l As before, we also include 
processed methods. 

We have not found methods of order eight. In fact, it is an open problem 
to determine if such a large system of polynomial equations admits solutions 
leading to more efficient methods than those collected in Table HI 

Runge Kutta Nystrom methods. As we have seen in section 13.71 meth- 
ods of this class may be considered as particular examples of composition (|36p . 
Nevertheless, their wide range of applicability to relevant physical problems has 
originated an exhaustive search of efficient schemes. Moreover, since in this case 
the associated vector fields F^ and F^ have different qualitative properties, 
methods with different features may be found in the literature. Thus, one may 
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Order 


3 


4 


6 


8 


3- [72] 


3-P31GSHES1E8! 

4-5-[5l] 


9- [33] 

10- ng 


27-? 




P:3,4-[II] 
P:2-7-[7] 


P:5-[II] 
P:5-10-[7] 


P:14? 



Table 5: Symmetric composition schemes of the form (|36p (appropriate when 
the ODE is split in two parts) and (|35p (composition of a method and its 
adjoint). At order eight, we have not found methods. They would require at 
least 27 stages or at least 14 stages for processed schemes. The notation is the 
same as in Tabled! 

find non-symmetric methods of the form 

AB = ip h = tp^ o tp^ h o • • • o tp^ h o ^ h 

BA = iph = <$ h o <f% h o ■ . • o $ h o ^ (83) 

where AS and are conjugate to each other, leading to the same perfor- 
mance. However, to take profit of the FSAL (First Same As Last) property, we 
can consider the following non equivalent compositions 

ABA ^ h = ^ +ih o V M o v W fc o.-o^o ^ (84) 

and 

= ^ fc = ^ ifc o v W fc o <$ h o ■ ■ ■ o o ^W fc . (85) 

The symmetric case (a s+ 2_i = aj,6 s+ i_j = 6j for the composition ABA and 
6 s+ 2-i = bi,a s+ i-i = a,i for the composition BAB) has also been proposed in 
this setting to get more efficient schemes. In this case, the minimum number 
of stages is, from Table [3] s r = L r — 1 (resp. k r = L* — 1), to get a method of 
order r (resp. effective order r). For non-symmetric compositions this minimum 
number can be obtained from Table [2j 

Highly efficient methods up to order six have been published. In Table [6] we 
collect the most representative within this class. We add S or N to distinguish 
symmetric from non-symmetric schemes and the subindex A B, ABA and BAB 
to denote compositions ([55]) . (JMJ) and (i85j) . respectively. Processed methods 
have also been included. 

To achieve order eight, the coefficients a^, bi in a non-processed scheme have 
to solve (in addition to consistency) a system of 16 nonlinear equations. A large 
number of solutions could exist, although, as far as we know, only one attempt 
to solve these equations has been reported [63] (the performance of such method 
was not clearly superior symmetric-symmetric methods). 

In [16] the authors have carried out a detailed analysis of the order conditions 
for symmetric compositions ABA and BAB. In this work 4th-order methods 
from 3 to 6 stages, and also 6th-order methods from 7 to 14 stages are analysed. 
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Order 


4 


5 


6 


8 


3S-IM1E01E81E8] 

4N AB -|56] 
4N BAB -H9] 

4-5Saba-M 
6Saba,bab-$M 


^ABA-m 
6^AB-M 


7S ABA -g3\&\ 
7Sbab-[33] 

8-15Saba,bab-[H] 

7,11S B ab-H51 


irs ABA -m 


P:2N AB -HI] 




P:^6Saba,bab-H31 
P:7S Byl B-[Il] 


P:9S ABA -m 
P:11S BA b-M 



Table 6: RKN splitting integrators. Since the role of the flows (p t and ip t 
is not interchangeable here, we distinguish symmetric S and non-symmetric N 
compositions with a subindex AB, ABA, BAB for the compositions (I83p . (|84p 
and ([85]) . As usual, processed methods are preceded by P. 



Order 


3 


4 


6 


8 


2N AB -(72] 


2Saba,bab-!1I1E6] 

4Saba,bab-[8Q] 
3,4Saba,bab-\M&\ 


4,5Saba,bab-[SI] 


HSABA,BAB-|67j 




P:1Sbab-M ED E3E] 
P:2S BA B-[5l] 


P^SylBA.BAB-lUJ 


P:4S A ba-[I3] 
P:5S Bj4b -[131[I4] 



Table 7: RKN splitting methods with modified potentials. Schemes are coded 
as in Table El 



The integrators selected perform extraordinarily well indeed. For instance, 
on the Henon-Heiles Hamiltonian (|19p the 4th-order 6-stage method is more 
accurate (at constant work) than leapfrog in a wide range of step sizes, whereas 
its global error is about 0.00175 times that of the classical 4th-order Runge- 
Kutta method. In consequence, its computational cost for a given error is about 
0.31. This has to be compared with the composition (|29p based on leapfrog, 
which have truncation errors about 10 times larger than the classical Runge- 
Kutta scheme. 

On the other hand, as we have seen in subsection l4.ll the particular structure 
of problem (|68|) allows one to use modified potentials in compositions (f83l) - 
(|85p . This is appealing when the evaluation of such modified potentials is not 
particularly costly. In such circumstances one may replace in (j83j) - (j85[ ) flows 
associated to hb{F^ with the corresponding to hCh(h, q, di, en, e^), as given in 
(1771) . The coefficients q, di, etc. can be used to solve some order conditions, so 
that methods with a reduced number of stages can be obtained. We emphasize 
that these schemes are of interest when the extra cost due to the modified 
potentials is moderate, as is the case in many problems arising in classical and 
quantum mechanics. In Table [7J we collect some relevant methods we have found 
in the literature, both processed and non-processed. 
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Order 


(n,2) 


(n,4) 


(n,5) 


l(2,2)S-p] 
n(2n,2)N A BA,BAB-\mM 


4(6,4)N BAB -P] 
5(8A)S A ba,bab-[53\ 




P: 1(32,2)- [87] 


P:3(7,6,4)S A ba-[12, 






P:2(6,4)Sab-[I2] 


P:3(7,6,5)S AB -[12] 




P:1(6A)S ABA -M 
P:n(n,4)S A BA-m\ 


P:2(7,6,5)S AB -H2] 



Table 8: Splitting methods for near-integrable systems. For processed methods 
we also include methods applicable when [[[-F' 6 ', F^], F^j, F^j = (second 
row) and schemes with modified flows (last two rows of processed methods). 



Methods for near-integrable systems. As we have seen in section 13.61 
splitting methods designed for equation (166j) have typically two relevant param- 
eters: h (the step size) and e (the size of the perturbation). In consequence, 
the dominant error in a given scheme depends on their relative size, and this 
depends usually on the particular problem considered (and sometimes even on 
the initial conditions). For this reason, a number of methods at different orders 
in both parameters h and e are found in the literature. We collect some of them 
in Table Here the notation is a bit clumsy: a method of order (n,4), say, 
means that the exact and the modified vector fields, i.e., hF and log(^(/i)) in 
(foTj) . differ in terms 0(eh n+1 + e 2 h 5 + •••), whereas in a method (7,6,4) this 
difference is 0(eh 8 + e 2 h 7 + e 3 /i 5 + •••). In both cases, the order of consistency 
in the limit h — > is four, but the last method incorporates more terms in the 
asymptotic expansion of the error. 

In [53] both families (ABA and BAB) of symmetric (2s, 2) schemes for 
s < 5 with positive coefficients are proposed which are about three times more 
accurate (at constant work) than leapfrog, whereas in [IB] a systematic study 
of (2s, 2) methods is carried out, obtaining new schemes up to s = 10 with 
positive coefficients. 

In some near- integrable problems, the identity [[[F®, F^], F^), F®] = 
still holds, where F^ is the vector field associated to /W, i = a,b, in (|66j) . This 
takes place, in particular, in Hamiltonian problems H = Hq + eH\ where Hq 
is quadratic in the kinetic energy and eH\ depends only on the coordinates 
(e.g. examples (|19p ~ (|22p can be split in this way, where Hq is the harmonic 
oscillator or the Kepler problem and H\ depends only on the coordinates). In 
consequence, the previous techniques used to obtain RKN methods still apply 
here, as well as the inclusion of flows of modified potentials in the composition. 

In Table [8] we separate, as usual, non-processed from processed schemes 
(preceded by P). In the later case we also include methods applicable when 
[[[-F 1 ^, FW], FM], FW] = (second row) and schemes with modified potentials 
(last two rows of processed methods). 
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6 Preserving properties and backward error analysis 



Much insight into the long-time behavior of splitting methods (including preser- 
vation of invariants and structures in the phase space) can be gained by applying 
backward error analysis techniques. We will summarize here some of the main 
issues involved and refer the reader to [5T] for a detailed treatment of the theory. 

When we analyzed in the Introduction the symplectic Euler scheme as ap- 
plied to the simple harmonic oscillator, we associated its good qualitative prop- 
erties with the fact that the numerical solution can be interpreted as the exact 
solution of a perturbed Hamiltonian system. This remarkable feature consti- 
tutes a simple illustration of the insight provided by backward error analysis 
(BEA) in this setting. More generally, suppose that we apply the splitting 
method (|36p to solve equation pip . Then the corresponding numerical solution 
at time t = h is given by 

x(h) = K(h)x = e b °+ lhB e a ° hA e b ° hB • • • e b2hB e aihA e blhB x , 

where the so-called stability matrix K(h) is given explicitly by 

1 \ / 1 a s h \ ( 1 ath W 1 
-b s+1 h 1 J \ 1 J "' I 1 J I -b x h 1 



K{h) = 
In this way, one gets 



K{h) 



K x {h) K 2 {h) 
K 3 (h) K 4 (h) 

where Ki(h), Ki{h) (respectively, K2(h), K 3 (h)) are even (resp. odd) func- 
tions and det-fC(Zi) = 1. As a matter of fact, any splitting method is uniquely 
determined by its stability matrix, so that the analysis can be carried out 
only with K(h) [9j. If in addition the splitting method is symmetric then 
K{h)~ l = K(—h) and we can write 



K{h) 



P(h) K 2 (h) 
K 3 (h) p(h) 



whevep(h) = \tx{K{h)) = \{KiQi) + K A {h)). It can be shown that the matrix 
K{h) is stable for a given h G R, i.e., K(h) n is bounded for all the iterations 
n, if and only if there exist real functions 4>(h),^(h) such that p(h) = cos((fr(h)) 
and K2{h) = —^(h) 2 K 3 (h). In that case 

/cos(WO) 7 (fc)sm(0(/O) \ / l(h)<l>(h)\ 

K ^ = {-^ cosim) J =exp l v -f} o J 

where, by consistency, ^'(0) = 1 and 7(0) = 1, whereas symmetry imposes 
<j){-h) = -<f>(h) and j(-h) = j(h). 

This result implies, in particular, that the numerical solution (q n ,Pn) at 
time t n = nh obtained by applying the splitting method to the linear system 
(fTTI) verifies 

cos(£ n d)) j(h) s'm(t n u;) \ / q 

-7(h)" 1 sin(t n a)) cos(t n a>) J \ Po 
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for values of h such that K{h) is stable. Here u = (fi(h)/h. Equivalently, 

Qn = q(tn), Pn = (</)(/l)7(/l)/^)~ 1 ^(in), 

where q(t) is the exact solution of 

d 2 , 

dT 2q + UJq = 

with initial condition q(0) = qo, q'(0) = (<j)(h)j(h)/h)po. In other words, the 
numerical solution provided by the splitting method is the exact solution of 
a harmonic oscillator with frequency u> ~ 1, i.e., of a system of equations 
satisfying the same geometric properties as the original system. The existence of 
such a backward error interpretation has direct implications for the qualitative 
behavior of the numerical solution, as well as for its global error. 

The main idea can be extended to an arbitrary non- linear ODE ([I]). Recall 
from Subsection 12.21 that each integrator i/j^ has associated a series = 
/ + h^i + h?"^2 + • • • of differential operators acting on smooth functions 
g 6 C°°(]R D ,]R), and its formal logarithm log(^(h)) is a series of vector fields 
(viewed as first order differential operators) log(^(h)) = hF\ + h 2 F 2 + • • • . 
For g £ C°°(]R D ,1R), the result of acting each on g is of the form Fk[g] = 
g'(x)fk(x), for a certain smooth map fk '■ — ► K . Now, consider the 
modified differential equation (defined as a formal series in powers of h) 

x' = f h (x) = f(x) + hf 2 (x) + h 2 h{x) + ■■■ (86) 

associated to the integrator vph. Then one has that x n = x(t n ), with t n = nh, 
which allows studying the long-time behaviour of the numerical integrator by 
analysing the solutions of the system (|86p viewed as a small perturbation of the 
original system ([TJ). This allows one to get important qualitative information 
about the numerical solution. In particular, 

• for symmetric methods, the modified differential equation only contains 
even powers of h; 

• for volume-preserving methods applied to a divergence-free dynamical 
system, the modified equation is also divergence-free; 

• for symplectic methods applied to a Hamiltonian system, the modified 
differential equation is (locally) Hamiltonian. 

In the particular case of a symplectic integration method, this means that 
there exist smooth functions Hj : M. 2d — > M for j = 2, 3, . . ., such that fj(x) = 
JVHj(x), where J is the canonical symplectic matrix. In consequence, there 
exists a modified Hamiltonian of the form 

H(q,p) = H(q,p) + hH 2 {q,p) + h 2 H 3 (q,p) + h*H A {q,p) + ■■■ (87) 

such that the modified differential equation is given by 

q' = V p H(q,p), p' = -V q H{q,p). 
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Of course, if the method has order r, say, then H \ = for i < r in (|87p . In 
other words, the modified Hamiltonian has the form H = H + h r H r+ i + • • • . 
In particular, for the Stormer-Verlet method (|l(jp applied to the Hamiltonian 
H(q,p) = T(p) + V(q), one has 



Apart from the linear case analyzed before, the series in ()86|) does not con- 
verge in general. To make this formalism rigorous, one has to give bounds on 
the coefficient functions fj(x) of the modified equation so as to determine an 
optimal truncation index and finally one has to estimate the difference between 
the numerical solution x n and the exact solution x{h) of the modified equation. 

These estimates constitute in fact the basis for rigorous statements about 
the long term behavior of the numerical solution. For instance, this theory 
allows one to proof rigorously that a symplectic numerical method of order 
r with constant step size h applied to a Hamiltonian system H verifies that 
H(x n ) = H(xq) + 0(h r ) for exponentially long time intervals |41j . 

On the other hand, since the modified differential equation of a numerical 
scheme depends explicitly on the step size used, one has a different modified 
equation each time the step size h is changed. This fact seems to be the reason 
of the poor long time behavior observed in practice when a symplectic scheme 
is implemented directly with a standard variable step-size strategy. 

7 Special methods for special problems 
7.1 Splitting methods for linear systems 

Suppose one is interested in solving numerically the differential equations aris- 
ing from the generalized harmonic oscillator with Hamiltonian function (|18|) . 
Although RKN methods with modified potentials can be always used for this 
purpose, we will see in the sequel that the particular structure of this system 
allows one to design specially tailored schemes which are orders of magnitude 
more efficient than other integrators frequently used in the literature. 

At this point, the reader could reasonably ask about the convenience of 
designing new numerical methods for the harmonic oscillator (|18p . It turns 
out, however, that efficient splitting methods for this system can be of great 
interest for the numerical treatment of partial differential equations appearing 
in quantum mechanics, optics and electrodynamics previously discretized in 
space. 

Suppose, in particular, that we have to solve numerically the time dependent 
Schrodinger equation (126j) with initial wave function tp(x,Q) = ipo{x). We can 
write ([26]) as 



H = H + h 2 ( — — 
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For simplicity, let us consider the one-dimensional problem and suppose that 
it is defined in a given interval x £ [xq,xn] (ip(xo,t) = ip(xN,t) = or it has 
periodic boundary conditions). A common procedure consists in taking first a 
discrete spatial representation of the wave function ijj(x,t): the interval is split 
in N parts of length Ax = (xjy — x )/N and the vector u = (u , . . . , un_i) t € 
is formed, with u n = if)(x n , t) and x n = xq + nAx, n = 0, 1, . . . , TV — 1. The 
partial differential equation (188j) is then replaced by the iV-dimensional linear 
ODE 

i— u (t) = H u(t), u(0) = u G C^, (89) 
at 

where H € ~R NxN represents the (in general Hermitian) matrix associated 
with the Hamiltonian [32]. The formal solution of equation (|89p is given by 
u(t) = e~** H uo, but to exponentiate this N x N complex and full matrix can 
be prohibitively expensive for large values of N, so in practice other methods 
are preferred. 

In general H = T + V, where V is a diagonal matrix associated with the 
potential energy V and T is a full matrix related to the kinetic energy T. Their 
action on the wave function vector is obtained as follows. The potential operator 
being local in this representation, one has (Vu) n = V(x n )u n and thus the prod- 
uct Vu requires to compute N complex multiplications. Since periodic bound- 
ary conditions are assumed, for the kinetic energy one has Tu = JF^Dy-^u, 
where T and J-~ l correspond to the forward and backward discrete Fourier 
transform, and Dy is local in the momentum representation (i.e., it is a diagonal 
matrix). The transformation T from the discrete coordinate representation to 
the discrete momentum representation (and back) is done via the fast Fourier 
transform (FFT) algorithm, requiring 0(N log N) operations. It is therefore 
possible to use the methods of subsection 12.11 with this splitting. 

There are other ways, however, of using splitting techniques in this context. 
To this end, notice that e~** H is not only unitary, but also symplectic with 
canonical coordinates q = Re(u) and momenta p = Im(u). Thus, equation 
(1891) is equivalent to [371 [38] 

|q = Hp, |p = "Hq, (90) 

where H q and H p require both a real-complex FFT and its inverse. In addi- 
tion, system (|90p can be seen as the classical evolution equations corresponding 
to the Hamiltonian function (|18p with M = N = H. Thus, efficient schemes 
for solving numerically the generalized harmonic oscillator can be applied di- 
rectly to this problem. Also the Maxwell equations (|28|) in an isotropic, lossless 
and source free medium, when they are previously discretized in space have a 
similar structure [70] . In consequence, numerical methods of this class are well 
adapted for their numerical treatment. 
Clearly, one may write 
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with the 2N x 2N matrices A and B given by 



The evolution operator corresponding to (|9ip is 

- f cos ( tH ) sin(tH) \ 
0(t) " V -sin(tH) cos(tH) J' (92) 

which is an orthogonal and symplectic 2Nx2N matrix. As before, its evaluation 
is computationally very expensive and thus some approximation is required. 
The usual procedure is to split the whole time interval into M steps of length 
h = t/M, so that 0(t) = [0(h)] , and then approximate 0(h) acting on the 
initial condition at each step. 
In this respect, observe that 

■*-(; ?)■ eB = ( -h i) 

and the cost of evaluating the action of e A and e B on z = (q, p) T is essentially 
the cost of computing the products H p and H q, respectively. It makes sense, 
then, to use splitting methods of the form (j36H . which in this context read 

O n (h) = e hb °+ lB e hasA ■ ■ ■ e hb2B e haiA e hblB . (93) 

Several methods with different orders have been constructed along these 
lines indeed [37\ [501 EH] • Of particular relevance are the schemes presented in 
[37j . since only s = r exponentials e haiA and e hb ' B are used to achieve order r for 
r = 4, 6, 8, 10 and 12. By contrast, in a general composition (|93|) the minimum 
number s of exponentials e haiA and e hbiB (or stages) required to attain order 
r = 8, 10 is k = 15,31, respectively [4T| |4"5]. 

Furthermore, one can use processing to reduce even more the number of 
exponentials. A different approach can also be followed, however: to take a 
number of stages larger than strictly necessary to solve all the order conditions 
to improve the efficiency and stability of the resulting schemes. The idea is to 
use the extra cost to reduce the size of the error terms, enlarge the stability 
interval and achieve therefore a higher efficiency but without raising the order. 

Kernels with up to 19, 32 and 38 stages have been proposed, and for each 
kernel the corresponding coefficients a,,6j have been determined according to 
two different criteria. The first set of solutions is taken so as to provide methods 
of order r = 10, 16 and 20. The second set of coefficients bring highly accurate 
second order methods with an enlarged domain of stability. A more detailed 
treatment can be found in (HJ [9]. 

7.2 Splitting methods for non-autonomous systems 

So far we have considered the problem of designing splitting methods for the 
numerical integration of autonomous differential equations ([T]). As we have 
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shown, there are a large number of schemes of different orders in the literature, 
and some of them are particularly efficient when the system possesses some 
additional structure, e.g., for the second-order differential equation y" = g(y) 
and the generalized harmonic oscillator (| 18[) . In this section we will review two 
different strategies to apply the splitting schemes when there is an explicit time 
dependency in the original problem. 

To fix ideas, let us assume that our system is non-autonomous and can be 
split as 

x' = f(x,t) = fW(x,t) + fW(x,t), x(0)=x . (94) 

The first, most obvious procedure consists in taking t as a new coordinate, 
so that (|94p is transformed into an equivalent autonomous equation to which 
standard splitting algorithms can be applied. More specifically, equation (f94j) 
is equivalent to the enlarged system 

d ( x \ j / M (^'«) 1 \ f lb \x,x a ) i 

*U«J I i J T o J 

V v ' V v ' 

/[i] ;p] 

with Xti,xt2 £ Note that if the systems 

y' = f [A] (y), y' = f [B] {y) 

with y = (x, xti, Xt2) are solvable, then a splitting method similar to (|36p can be 
used, since Xti is constant when integrating the first equation and xt2 is constant 
when solving the second one. This, in fact, can be considered as a generalization 
of the procedure proposed in [73] for time-dependent and separable Hamiltonian 
systems, and is of interest if the time-dependent part in and is cheap to 
compute. Otherwise the overall algorithm may be computationally costly, since 
these functions have to be evaluated s times (the number of stages in (f36j) ) per 
time step. 

Another disadvantage of this simple procedure is the following. Suppose 
that, when the time is frozen, the function / in (|94p has a special structure 
which allows to apply highly efficient splitting schemes. If now t is a variable, 
with ([95]) this time dependency is eliminated but the structure of the equation 
might be modified so that one is bound to resort to more general and less 
efficient integrators. This issue has been analyzed in detail in [5]. 

A second procedure which avoids the difficulties exhibited by the previous 
example consists in approximating the exact solution of (|94p or equivalently the 
flow <fh by the composition 

llr] [Bs+i] [As] [B a ] [B 2 ] [At] [Bi] / nft x 

%,h = < Ph ol Ph 0i Ph J , (96) 

where the maps <p\?*\ are ^ ne exac t 1-flows corresponding to the time- 

independent differential equations 

x' = Ai(x), x' = Bi(x), i = l,2,... (97) 
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(95) 



respectively, with 

k k 
Ai(x) ^h^Pijf^i^-rj), = hJ2 a vf [b] ( x > T i)- ( 98 ) 

3=1 3=1 

Here Tj = to + c jh an d the (real) constants Cj, pij, Oij are chosen such that 
= V ; i r l + C(/i r+1 ). Furthermore, the new schemes, when applied to (|94p with 
the time frozen, reproduce the standard splitting (|36p . This is accomplished by 
ensuring that Ylj Pij = a i anc ^ Ylj a ij = &t- The Cj coefficients, on the other 
hand, are typically chosen as the nodes of a symmetric quadrature rule of order 
at least r. In particular, if a Gauss-Legendre quadrature rule is adopted, with 
k evaluations of (x, tj) and (x, tj) a method of order r = 2k can be built 
(taking s sufficiently large). 

Once the quadrature nodes Tj and the number of stages s are fixed, there still 

M 

remains to obtain the coefficients pij, aij such that tp s h has the desired order. 
This is done by requiring that the composition (I96j) match the solution of (|94j) 
as given by the Magnus expansion [TTJ . The task is made easier by noticing that 
the order conditions to be satisfied by and are identical both for linear 
and nonlinear vector fields. Thus, the problem for the linear case is solved 
first and then one generalizes the treatment to arbitrary nonlinear separable 
problems. 

The integrators of order four and six constructed along these lines in [5] are 
generally more efficient than standard splitting methods applied to the enlarged 
system (|95|) . 

8 Numerical examples with selected methods 

This section intends to illustrate the relative performance between different 
splitting methods, and occasionally we compare with other standard methods. 
We consider first a relatively simple problem where most of the methods pre- 
viously mentioned can be used, showing their good features. We show the 
interest of the high order methods when accurate results are desired and the 
improvement which can be achieved when choosing a method from the most 
appropriate family of methods for each problem. Next, we consider a problem 
which, due to its very particular structure, allows to build tailored methods 
whose performance is much superior to other splitting methods. 

8.1 The perturbed Kepler problem 

As a first example, we take the perturbed Kepler problem with Hamiltonian 
(HQ} 

where r = \J q\ + an d the additional parameter a has been introduced for 
convenience. This Hamiltonian describes in first approximation the dynamics 
of a satellite moving into the gravitational field produced by a slightly oblate 
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spheric planet. The motion takes place in a plane containing the symmetry axis 
of the planet when a = 1, whereas a = corresponds to a plane perpendicular 
to that axis [60J. 

This simple (but non trivial) example constitutes in fact an excellent test 
bench for most of the methods of this paper. Notice that the system is separable 
into kinetic and potential parts, and we can use, for instance, the symmetric 
second order method (|10j) which allows us to get higher order methods by 
composition, as given in (|62p . On the other hand, since the system is separable 
into two solvable parts, then we can also use methods from Table El which 
should show better performances than methods of the same order considered 
from the previous family of methods. In addition, the kinetic energy is quadratic 
in momenta, so that RKN methods from Table [6] can be used, and one expect 
a further improvement. Finally, observe that one may split the system as 

H = Ho + eH h (100) 

where Hq corresponds to the Kepler problem, which is exactly solvable. The 
Keplerian part of the Hamiltonian can be solved in action- angle coordinates, 
where two changes of variables are needed. Alternatively, if desired, Hq can be 
integrated in cartesian coordinates using the / and g Gauss functions, but then 
a nonlinear equation must be solved with an iterative scheme [31]. In any case, 
if e <C 1, methods from Table [8] can be used which should be superior to all 
previous methods in the limit s — > 0. 

We must also mention that the performance of all methods previously men- 
tioned can be further improved by using the processing technique, and even 
additional improvements can be achieved if modified potentials are considered. 

We take e = 0.001, which approximately corresponds to a satellite moving 
under the influence of the Earth [46] and initial conditions q\ = 1 — e, qi = 0, 
pi = 0, pi = a/(1 + e)/(l — e), with e = 0.2 (this would be the eccentricity 
for the unperturbed Kepler problem). In general, no closed orbits are present 
and a precession is observed. Notice that for the Hamiltonian (|99p the strength 
of the perturbation depends obviously of the value of e, but also on the initial 
conditions. We take a = 1 and determine numerically the trajectory for up 
to the final time tf = 500 • 2ir (the exact solution is accurately approximated 
using a high order method with a very small time-step, and this computation 
was repeated with different time steps and methods to assure the accuracy is 
reached up to round off). 

To compare the performance of different methods it is usual to consider 
efficiency curves. We measure the average error in energy computed at times 
tk = k ■ 2tt for k = 401, 402, . . . , 500 and this is repeated several times for each 
method and using different time steps (changing the computational cost for the 
numerical integration) . 

In the first numerical test, we compare the relative performance between 
different symmetric-symmetric methods collected in Table HI We choose as 
the basic method the symmetric second order composition (jlOj) to build higher 
order methods with the composition (|62j) . As mentioned, in general, the per- 
formance of the methods of the same order increase with the number of stages 
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for the methods in Table [H This is illustrated in Figure. H] where we show the 
performance of two 4th-order methods with three stages (given by (|29j) ) and 
five stages (given by (|82j) ). The results show that for this problem the five stage 
method is more accurate for all computational costs considered. A similar fea- 
ture is observed for the other methods at higher orders (with the exception of 
the 24-stage 8th-order method which was obtained in a different way and the 
21-stage methods shows a better performance). We choose the best method 
from Table H] at each order (including the well known three-stages 4th-order 
method as a reference) where we denote by SS s r the corresponding method of 
order r using a s-stage composition: 



SSi2: The 2nd-order method (|10|) which has the highest possible stability 
among splitting methods. 



• SS34 The well known 3-stage 4th-order method <\29h . 

• SS 5 4 The 5-stage 4th-order method dHSJ) [78]. 

• SS136, SS218, SS35IO: The composition from Table [4] and whose coeffi- 
cients are given in |76j . 

The results are shown in Fig. HI where we clearly observe that the high order 
methods have better performance when high accuracy is desired. 
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Figure 4: Average error in energy versus number of force evaluations in a dou- 
ble logarithmic scale for the numerical integration of the Hamiltonian system 
(|99p . It is shown performance of the most efficient non-processed symmetric- 
symmetric methods from Table 0) 

The following numerical experiment intends to illustrate the interest of the 
methods designed for problems with some particular structure. For simplicity, 
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in this numerical test we only consider fourth-order methods from different 
families of methods which can be used on this problem, in order to observe the 
benefit of tuned methods for problems with particular structures. The following 
methods are considered in addition to SS34 and SSs4: 

• Sg4: The symmetric 6-stage 4th-order method for separable problems |16j 
from Table [5j 

• RKNg4: The symmetric 6Sbab 4th-order method for Nystrom problems 
[16] from Table El 

• NI(8,4): The 5-stage fourth-order method 5(8,4)Sbab given in [53] from 
Table El 

• RK44: The standard 4-stage 4th-order non-symplectic Runge-Kutta meth- 
ods, used reference method. 

Figure [5] shows in double logarithmic scale the results obtained. In the left 
panel we show the average error in energy versus the number of force evaluations 
and in the right panel we repeated the same experiment, but we measured the 
average error in position (computed at the same instants). For the method 
NI(8,4) this counting of the computational cost is not an appropriate measure. 
Its computational cost strongly depends on each particular problem since the 
evolution of Hq has to be computed exactly (or very accurately) . For simplicity, 
in our experiments, we have considered that one stage of NI(8,4) is twice as 
expensive as one evaluation of the force. We have also included as a reference 
the curve obtained in Fig. H]by SS35IO. 

Observe that in the first case the results will be largely independent of tf 
because the average error in energy does not increase secularly for symplectic 
integrators. For comparison, we have also included the results obtained by 
the standard 4-stage fourth-order Runge-Kutta method whose error in energy 
grows linearly and the error in positions quadratically. 

Even more accurate results could be obtained as follows. As mentioned, for 
this particular problem, modified potentials could be used and this can be done 
at a very low computational cost. Then, methods from Tables \7\ and [8] can be 
used. For instance, for the split (jlOOp we can apply methods which incorporate 
modified perturbations exp(eC^(6, c)) into the algorithm. Then the following 
map has to be evaluated: 



where A = (3/2)(a(3g 2 - 2g 2 ) - r 2 ), B = (3/2)(a5g 2 - r 2 ), C = 9(2r 4 + 
3ar 2 {ql-Aql) + a 2 {lSq\ + qlql-2q%)) and D = 9(2r 4 - 15ar 2 g 2 + 5a 2 g 2 (5g 2 + 



2g|)). Notice that the increment in the computational cost with respect to the 




(101) 



evaluation of e hebF (which corresponds to c = 0) is only due to a few very 



44 



SS„10 



-12 



-6 



4.4 



4.6 



4.8 5 5.2 5.4 

LOG(N. EVALUATIONS) 



5.6 



4.4 



4.6 



4.8 5 5.2 '. 

LOG(N. EVALUATIONS) 



5.4 



5.6 



Figure 5: Average error in energy (left panel) and position (right panel) versus 
number of force evaluations in a double logarithmic scale for the numerical 
integration of the Hamiltonian system (I99j) . The performance of different 4th- 
order methods from Tables 0H8] is shown. As a reference, we also shown the 
results obtained by the standard 4-stage fourth-order Runge-Kutta method. 

simple additional operations. For this particular example, the evaluation of the 
modified perturbation e eC h(b,c) - 1S ^out a 10 — 20% more expensive than e hbF ^ . 

As a result, more elaborated and efficient methods are obtained by consid- 
ering the schemes (n, 4) and (n, 5) from Table [8] with processing (which only 
requires a few more code lines to program) for RKN problems and which incor- 
porate modified potentials (see [T2]). 

8.2 The Schrodinger equation 

As a second example we consider the one-dimensional time-dependent Schro- 
dinger equation (|26p with the Morse potential V(x) = D (1 — e~ ax ) 2 . We fix 
the parameters to the following values in atomic units (a.u.): fi = 1745 a.u., 
D = 0.2251 a.u. and a = 1.1741 a.u., which are frequently used for modelling 
the HF molecule. As initial conditions we take the Gaussian wave function 
V>(x,i) = /oexp ( — (3{x — x) 2 ), with (3 = \fkjl/2, k = 2Da 2 , x = —0.1 and p 
is a normalizing constant. Assuming that the system is defined in the interval 
x £ [—0.8, 4.32], we split it into d = 128 parts of length Ax = 0.04, take periodic 
boundary conditions and integrate along the interval t S [0,20 • 27r/u;o] with 
f^o = d\JlD I p. (see [8] for more details on the implementation of the splitting 
methods to this particular problem). 

Figure [6] shows the error in the Euclidean norm of the vector solution at the 
end of the integration versus the number of FFT calls in double logarithmic 
scale. The integrations are done starting from a sufficiently small time step and 
repeating the computation by slightly increasing the time step until an overflow 
occurs, which we identify with the stability limit. We present the results for 
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the following methods (in addition to the previous ones SSi2, SS34 and SS218): 

• RKNn6: The 11-stage 6th-order method 11Sr4r-|16] from Table 

• GM12I2: The 12-stage 12th-order method from [37] tailored for linear 
problems with this particular structure. 

• P382: The 38-stage second order processed method with coefficients given 
in [8] tailored for linear problems with this particular structure (only the 
computational cost required to evaluate the kernel has been taken into 
consideration) . 

• T r r, r = 8, 12: r-stage rth-order Taylor methods obtained by truncating 
the exponential up to order r. 
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Figure 6: Error of the vector solution for the Schrodinger equation versus the 
number of FFT calls in a log- log scale for the symmetric-symmetric composition 
methods: SS^n, for methods of order n using /e-stage compositions; RKNn6 is 
an 11-stage 6th-order methods from [16] for Nystrom problems; the 12-stage 
12th-order method, GM12I2, from [37J; and P382, a 38-stage second order pro- 
cessed method. 

From the figure we observe that, for this numerical experiment, standard 
Taylor methods outperform to general symmetric-symmetric splitting methods 
and are also more accurate than the RKN method. This is because this problem 
has a very particular structure and these splitting methods are not optimized 
for them. However, the schemes GM12I2 and P382 are built for linear problems 
with this structure and their superiority is clearly manifest. It is important to 
remember that Taylor methods are non-geometric integrators. The numerical 
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experiments are carried out for a relatively short time, and the relative perfor- 
mances of the Taylor methods deteriorates with respect to splitting methods 
for longer integrations. 

9 Conclusions and outlook 

Splitting methods are a flexible and powerful way to solve numerically the ini- 
tial value problem defined by ([1]) when / can be decomposed into two or more 
parts and each of them is simpler to integrate than the original problem. This 
is especially true when the exact flow possesses some structural features which 
seems natural to reproduce at the discrete level, as happens, for instance, in 
Hamiltonian, Poisson, volume-preserving or time-reversible dynamical systems. 
They are explicit, usually simple to apply in practice and constitute an impor- 
tant class of geometric numerical integrators. Closely connected with splitting 
schemes are composition methods. In this case, the idea is to construct numeri- 
cal integrators of arbitrarily high order by composing one or more basic schemes 
of low order with appropriately chosen coefficients. The resulting method in- 
herits the relevant properties that the basic integrator shares with the exact 
solution, provided these properties are preserved by composition. 

In this paper we have reviewed some of the main features of splitting and 
composition methods in the numerical integration of ordinary differential equa- 
tions. We have presented a novel approach to get the order conditions of this 
class of schemes based on Lyndon words and we have seen how these order 
conditions particularize when coping with special classes of dynamical systems 
(near-integrable systems and second-order differential equations of the form 
y" = 9(y))- If turns out that the number of equations to be solved increases 
dramatically with the order considered, as so does the complexity of the prob- 
lem of finding efficient high order methods. One way to circumvent (up to a 
certain point) this difficulty consists in applying the processing technique, since 
then it is possible to design algorithms with fewer evaluations per time step. 
In this sense, one could say that the use of processing is perhaps the most 
economical path to achieve high order. 

Since splitting methods are widely applied in many areas of science, it is 
not surprising that a great number of different schemes are available in the 
mathematical, physical and chemical literature. We have collected here some 
of the most representative integrators, classified according to the particular 
structure of the differential equations, the number of stages and the order of 
consistency, citing in each case the actual reference where the method has been 
first proposed. 

The good qualitative behavior exhibited by splitting methods (including 
preservation of invariants and structures in phase space), as well as their fa- 
vorable error propagation in long-time integrations can be accounted for by 
applying the theory of backward error analysis. Loosely speaking, the observed 
performance is related with the fact that the numerical solution provided by the 
splitting method is the exact solution of a differential equation with the same 
geometric properties as the original system. This interpretation constitutes in 
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addition the basis for rigorous estimates on the numerical solution. 

In contrast with standard integration methods (Runge-Kutta, multistep), 
whose efficiency is essentially independent of the particular differential equation 
considered, splitting schemes can be designed to incorporate in their formulation 
some of the most relevant properties of the original system. This feature has 
to be taken into account when comparing the efficiency of splitting methods 
with respect to other general purpose integrators. In this sense, Figures [4] and 
[5] are quite illustrative. For this particular problem, specially adapted 4th- 
order splitting schemes are up to 6 orders of magnitude more accurate with 
the same computational cost than the well known Runge-Kutta method. They 
even outperform other standard higher order composition integrators for a wide 
range of values of the step size h. 

As an additional evidence of the extraordinary flexibility of splitting meth- 
ods, we have considered the problem of designing specially tailored schemes for 
the numerical integration of the generalized harmonic oscillator (|18[) . It turns 
out that several partial differential equations appearing in quantum mechanics, 
optics and electrodynamics give rise, once discretized in space, to this system 
with different matrices M and N. The particular structure of this dynamical 
system can be exploited to build an optimized processed second order method 
involving a large number of stages that nevertheless is far more efficient than 
other integrators. 

There are other issues in connection with splitting and composition methods 
that we have not tackled here, however, and that are also important in this 
context. Among them we can mention the following. 

• As was remarked in the introduction, no general rule is provided here to 
split any given function / in the differential equation (pQ). It turns out 
that, for / within a certain class of ODEs, this can be done systematically, 
whereas for other functions one has to proceed on a case by case basis. 
Sometimes, several splittings are possible, and the different schemes built 
from them lead to the preservation of distinctive geometric properties. It 
makes sense, then, to classify the ODEs and their corresponding integra- 
tion methods into different categories. This aspect has been analyzed in 
|57| . Moreover, in many physical problems there are several geometric 
properties that are conserved simultaneously along the evolution and it 
is not clear at all how to design methods preserving all of them. In that 
case, which one is the most relevant from a numerical point of view? 

• In this paper we have only considered the initial value problem defined 
by eq. dH) and integration methods with constant step size h. Backward 
error analysis provides an argument why this has to be the case in geo- 
metric numerical integration: the modified equation corresponding to the 
numerical method depends explicitly on h, so that if h is changed so does 
the modified equation and no preservation of geometric properties is guar- 
anteed. There are problems, however, where the use of an adaptive step 
size is mandatory, for instance in configurations of the iV-body problem 
allowing close encounters. In this case one may apply splitting methods 
with variable step size by using some specifically designed transformations 
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involving the time variable, in such a way that in the new variables the 
resulting time step is constant (see, e.g., [3]). 

• As we have shown in section 13.51 the presence of negative coefficients in 
splitting methods of order higher than two is unavoidable. This is not 
a problem when the flow of the differential equation evolves in a group 
(such as in the Hamiltonian case), but may be unacceptable when the 
ODE originates from a partial differential equation that is ill-posed for 
negative time progression. Several alternatives have been proposed in the 
literature, mainly by considering, when possible, modified potentials [3], 
as noted in section 14.11 One should observe, however, that the analysis 
done in section 13.51 does not preclude the existence of complex coefficients 
with positive real parts. As a matter of fact, splitting methods with com- 
plex coefficients have been developed and tested for problems in which the 
Hamiltonian is split into kinetic and potential energy terms [23J, for the 
time-dependent Schrodinger equation [2], for generic parabolic equations 
[23] and also in the more abstract setting of evolution PDEs in analytic 
semigroups [42]. 

• An important characteristic of any numerical integration method is sta- 
bility. Roughly speaking, the numerical solution provided by a stable 
numerical integrator does not tend to infinity when the exact solution is 
bounded. Although important, this feature has received considerably less 
attention in the specific case of splitting methods. To test the (linear) 
stability of the method (j36|) . instead of the linear equation y' = ay as in 
the usual stability analysis for ODE integrators, one considers the har- 
monic oscillator y" + \ 2 y = 0, A > 0, as a model problem with a splitting 
of the form (jlip . The idea is to find the time steps for which all numerical 
solutions remain bounded. The integrator (I36j) typically will be unstable 
for \hX\ > x*, where the parameter x* determines the stability threshold 
of the numerical scheme. In particular, for the leapfrog method one has 
x* = 2. Although the stability threshold imposes restrictions on the step 
size, in the process of building high order schemes, linear stability is not 
usually taken into account, ending sometimes with methods possessing 
such a small relative stability threshold that they are useless in practice. 
In this way, constructing high order splitting methods with relatively large 
linear stability intervals and highly accurate is of great interest. This has 
been achieved in reference [9] for linear systems, but remains an open 
problem in general. 

• In section [5] we have mentioned an optimization criterion to choose the free 
parameters in splitting and composition methods, which consist in mini- 
mizing the Euclidean norm of the coefficients that constitute the leading 
error term of the method. It is clear, however, that minimizing the lead- 
ing error term does not guarantee that the method thus obtained is the 
most efficient: it might occur that the influence of the subsequent error 
terms is the decisive factor in the performance of the scheme. In this 
sense, it would be extremely interesting to have estimates on all the error 
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terms in the asymptotic expansion of the modified equation and get the 
coefficients of the method that minimize these estimates. 

• The numerical analysis of second-order differential equations with oscilla- 
tory solutions has aroused much interest during the past few years. The 
typical test problem in this setting is the equation q" + Q 2 q = f(q), where 
is a symmetric and positive definite matrix. Here the aim is to de- 
sign new methods which improve in accuracy and stability the standard 
Stormer-Verlet integrator. We refer the reader to [41] and references 
therein for a comprehensive study of this problem. 

• Although only ODEs have been considered here, splitting methods have 
been also applied with success to stochastic differential equations (SDEs). 
Here the aim is, as in the deterministic case, to design integration meth- 
ods which automatically incorporate conservation properties the SDE pos- 
sesses [6l] . 
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